Linking community science data with population dynamics under stochastic process in the context of ecological release
Orlando Acevedo-Charry, José Miguel Ponciano, Scott K. Robinson, & Miguel A. Acevedo - Natya Hans is supporting HiperGator opperations and will be more involved for publications
2024-04-17
Could eBird data track population dynamics in time-series in the context of ecological release?
Ecological release: Population level niche expansion and/or shift when antagonistic interspecific interaction is removed or reduced (Herrmann et al., TREE 2021). Hypothetical ecological release could be explored as a function of species richness. Lower competition in islands expand niche (e.g., abundance), while more habitats and more species increase level of specialization in the Mainland (Sherry et al., Auk 2020).
Using the Caribbean meta-archipelago, we filtered the no-aquatic birds reported in eBird, between Julian day 125 and 150 from 2012-2013, between 6 and 11 hours in the morning. We excluded incomplete list and retained checklists of only “traveling” or “stationary” protocol, with a minimum effort of less than 2 km, more than 10 minutes but less than 2 hours (10:120 min) and less than 10 observers. We assumed as counting samples the mean of counts of each species reported by the observers during the 25 days filtered per year; if the mean is very constant among independent observers, we assume a higher probability of detection, while higher variation of the counts might reflect a difference in detection probability. These data constitutes the primary time-series of species for different islands.
To test a model of stochasticity process, we selected a widespread and easy to detect species in the meta-archipelago, which are distributed from the Mainland to different islands in the Caribbean: Bananaquit Coereba flaveola. Theory predicts that species in the mainland have lower density than islands, and that density could increase as species richness and potentially interspecific competition decreases.
Filtering eBird data in meta-archipelagos
The code in HiperGator was (Do not run):
########################################
###~ eBird data - Island filtering ~####
########################################
#Packages
library(auk); #eBird data filters
library(tidyverse); #manipulate data
library(dggridR) #Spatiotemporal Subsampling
library(iNEXT) #rarefaction
#~ eBird data ~####
#Data downloaded from eBird, up to August 2023, saved in the folder of the HPG
#Define the range
CaribbeanRange <- c(-91,8,-59,27)
IndoMalayRange <- c(90,-20,180,20)
#select the columns to extract (based on x$col_idx$name)
colsE <- c("observer_id", "sampling_event_identifier",
"group identifier",
"common_name", "scientific_name",
"observation_count",
"country", "state_code", "locality_id", "latitude", "longitude",
"protocol_type", "all_species_reported",
"observation_date",
"time_observations_started",
"duration_minutes", "effort_distance_km",
"number_observers")
#For Caribbean
f_ebdCa <- "ebd_IslandsCa.txt" #Temporal file to save the filtering eBird data (records)
f_sedCa <- "sed_IslandsCa.txt" #Temporal file to save the filtering sampling event data
ebd_filtCarib<-auk_ebd("eBirddata/ebd_relAug-2023.txt",
file_sampling = "eBirddata/ebd_sampling_relAug-2023.txt") %>%
auk_bbox(CaribbeanRange) %>% #W, S, E, N
auk_year(c(2012:2023)) %>%
auk_protocol(c("Traveling", "Stationary")) %>%
auk_distance(distance = c(0,5)) %>%
auk_duration(duration = c(0,300))%>%
auk_complete() %>%
auk_filter(f_ebdCa, f_sedCa, overwrite=T, keep = colsE)
f_ebdIM <- "eBirddata/ebd_IslandsIM.txt" #Temporal file to save the filtering eBird data (records)
f_sedIM <- "eBirddata/sed_IslandsIM.txt" #Temporal file to save the filtering sampling event data
ebd_filtIndoP<-auk_ebd("eBirddata/ebd_relAug-2023.txt",
file_sampling = "eBirddata/ebd_sampling_relAug-2023.txt") %>%
auk_bbox(IndoMalayRange) %>% #W, S, E, N
auk_year(c(2012:2023)) %>%
auk_protocol(c("Traveling", "Stationary")) %>%
auk_distance(distance = c(0,5)) %>%
auk_duration(duration = c(0,300))%>%
auk_complete() %>%
auk_filter(f_ebdIM, f_sedIM, overwrite=T, keep = colsE)
#and with read_ebd I apply another filter to do not repeat records from groups
sed_Carib <- read_sampling(f_sedCa)
ebd_Carib <- read_ebd(f_ebdCa)
sed_IndoM <- read_sampling(f_sedIM)
ebd_IndoM <- read_ebd(f_ebdIM)
# Function to convert time observation to hours since midnight
time_to_decimal <- function(x) {
x <- hms(x, quiet = TRUE)
hour(x) + minute(x) / 60 + second(x) / 3600
}
#Caribbean data
Carib_eff <- ebd_Carib %>%
mutate(
# I don't have here count in X to convert to NA
observation_count = as.integer(observation_count),
# effort_distance_km to 0 for non-travelling counts
effort_distance_km = if_else(protocol_type == "Stationary",
0, effort_distance_km),
# convert time to decimal hours since midnight
time_observations_started = time_to_decimal(time_observations_started),
hour_sampling = round(time_observations_started, 0),
# split date into year, month, week, and day of year
year = year(observation_date),
month = month(observation_date),
week = week(observation_date),
day_of_year = yday(observation_date)) %>%
filter(number_observers <= 10, #Only list with less than 10 observers
effort_distance_km <= 5, #be sure of distance effort
duration_minutes %in% (0:300), #be sure of duration effort
!is.na(observation_count)) #only records with abundance estimation
Checklists_with_more_than_fourSppC <- Carib_eff %>%
select(sampling_event_identifier, scientific_name) %>%
group_by(sampling_event_identifier) %>%
summarise(N_species = n()) %>%
filter(N_species >= 4)
#352797 lists
Caribbean_cooc <- Carib_eff[Carib_eff$sampling_event_identifier %in% Checklists_with_more_than_fourSppC$sampling_event_identifier, ]
Carib_sed <- sed_Carib %>%
mutate(
# effort_distance_km to 0 for non-travelling counts
effort_distance_km = if_else(protocol_type == "Stationary",
0, effort_distance_km),
# convert time to decimal hours since midnight
time_observations_started = time_to_decimal(time_observations_started),
hour_sampling = round(time_observations_started, 0),
# split date into year, month, week, and day of year
year = year(observation_date),
month = month(observation_date),
week = week(observation_date),
day_of_year = yday(observation_date)) %>%
filter(number_observers <= 10, #Only list with less than 10 observers
effort_distance_km <= 5, #be sure of distance effort
duration_minutes %in% (0:300)) #be sure of duration effort
#And now with the Oriental-Indo_Malayan-Papuan-Melanesian
IndoM_eff <- ebd_IndoM %>%
mutate(
# I don't have here count in X to convert to NA
observation_count = as.integer(observation_count),
# effort_distance_km to 0 for non-travelling counts
effort_distance_km = if_else(protocol_type == "Stationary",
0, effort_distance_km),
# convert time to decimal hours since midnight
time_observations_started = time_to_decimal(time_observations_started),
hour_sampling = round(time_observations_started, 0),
# split date into year, month, week, and day of year
year = year(observation_date),
month = month(observation_date),
week = week(observation_date),
day_of_year = yday(observation_date)) %>%
filter(number_observers <= 10, #Only list with less than 10 observers
effort_distance_km <= 5, #be sure of distance effort
duration_minutes %in% (0:300), #be sure of duration effort
!is.na(observation_count)) #remove records without abundance estimation
Checklists_with_more_than_fourSppIP <- IndoM_eff %>%
select(sampling_event_identifier, scientific_name) %>%
group_by(sampling_event_identifier) %>%
summarise(N_species = n()) %>%
filter(N_species >= 4)
IndoM_eff_cooc <- IndoM_eff[IndoM_eff$sampling_event_identifier %in% Checklists_with_more_than_fourSppIP$sampling_event_identifier, ]
IndoM_sed <- sed_IndoM %>%
mutate(
# effort_distance_km to 0 for non-travelling counts
effort_distance_km = if_else(protocol_type == "Stationary",
0, effort_distance_km),
# convert time to decimal hours since midnight
time_observations_started = time_to_decimal(time_observations_started),
hour_sampling = round(time_observations_started, 0),
# split date into year, month, week, and day of year
year = year(observation_date),
month = month(observation_date),
week = week(observation_date),
day_of_year = yday(observation_date)) %>%
filter(number_observers <= 10, #Only list with less than 10 observers
effort_distance_km <= 5, #be sure of distance effort
duration_minutes %in% (0:300)) #be sure of duration effort
#Only those land birds ####
#AVONET filtered by No aquatics
filterAvonet <- read.csv("filterAVONET/AVONET_filter_NoAquatics.csv")
#to combine the data, we have to adjust the names
names(filterAvonet)[2] <- "scientific_name"
#filter in eBird data those species identify in AVONET as land birds
ebird.carib = Caribbean_cooc[Caribbean_cooc$scientific_name %in% filterAvonet$scientific_name, ]
ebird.indop = IndoM_eff_cooc[IndoM_eff_cooc$scientific_name %in% filterAvonet$scientific_name, ]
##Spatial Grids ####
set.seed(123)
# Spatiotemporal Subsampling in diameters of ~1 km (area of ~100 km^2)
dggs_island <- dgconstruct(spacing = 11) #spacing 1 correspond to Characteristic Length Scale, or diameter of spherical cell
#add a new variable that identify cell
Caribbean_SS <- ebird.carib %>%
mutate(cell = dgGEO_to_SEQNUM(dggs_island, #id for cells
longitude, latitude)$seqnum) %>%
group_by(cell, scientific_name) %>%
mutate(mu_count = mean(observation_count))
IndoPacific_SS <- ebird.indop %>%
filter(observation_count < 7000) %>% #there is an outlying count for _Eurylaimus ochromalus_ https://ebird.org/checklist/S112199250
mutate(cell = dgGEO_to_SEQNUM(dggs_island, #id for cells
longitude, latitude)$seqnum) %>%
group_by(cell, scientific_name) %>%
mutate(mu_count = mean(observation_count))
#Get the number of species in each cell
RichnessCellCaribbean <- Caribbean_SS %>%
group_by(cell, scientific_name) %>%
summarise(count=n()) %>%
group_by(cell) %>%
summarise(SpRichness=n()) %>%
filter(SpRichness >= 4)
#Exclude cells in SS with less than 4 species
Caribbean_SS = Caribbean_SS[Caribbean_SS$cell %in% RichnessCellCaribbean$cell, ]
#Community assemblages
Caribbean_Cell_Assembly <- Caribbean_SS %>%
select(cell, scientific_name, mu_count) %>%
summarise(mu_count = mean(mu_count)) %>%
pivot_wider(names_from = cell, values_from = mu_count, values_fill = 0)
#Get the grid cell boundaries for cells which had quakes
gridCaribbean <- dgcellstogrid(dggs_island,RichnessCellCaribbean$cell)
#Update the grid cells' properties to include the number of lists in each cell
gridCaribbean <- merge(gridCaribbean,RichnessCellCaribbean,by.x="seqnum",by.y="cell")
# Handle cells that cross 180 degrees
wrapped_gridCaribbean = st_wrap_dateline(gridCaribbean,
options = c("WRAPDATELINE=YES","DATELINEOFFSET=180"), quiet = TRUE)
#Backup
saveRDS(wrapped_gridCaribbean, "Completeness_data/wrapped_gridCaribbean.rds")
saveRDS(Caribbean_SS, "Completeness_data/Caribbean_SS_140324.rds")
saveRDS(Caribbean_Cell_Assembly, "Completeness_data/Caribbean_Cell_Assembly.rds")
#Get the number of records in each cell
RichnessCellIndoPacific <- IndoPacific_SS %>%
group_by(cell, scientific_name) %>%
summarise(count=n()) %>%
group_by(cell) %>%
summarise(SpRichness=n()) %>%
filter(SpRichness >= 4)
#Exclude cells in SS with less than 4 species
IndoPacific_SS = IndoPacific_SS[IndoPacific_SS$cell %in% RichnessCellIndoPacific$cell, ]
#Community assemblages
IndoPacific_Cell_Assembly <- IndoPacific_SS %>%
select(cell, scientific_name, mu_count) %>%
summarise(mu_count = mean(mu_count)) %>%
pivot_wider(names_from = cell, values_from = mu_count, values_fill = 0)
#Get the grid cell boundaries for cells which had quakes
gridIndoPacific <- dgcellstogrid(dggs_island,RichnessCellIndoPacific$cell)
#Update the grid cells' properties to include the number of lists in each cell
gridIndoPacific <- merge(gridIndoPacific,RichnessCellIndoPacific,by.x="seqnum",by.y="cell")
# Handle cells that cross 180 degrees (and some weird in the map)
wrapped_gridIndoPacific = st_wrap_dateline(gridIndoPacific,
options = c("WRAPDATELINE=YES","DATELINEOFFSET=180"), quiet = TRUE)
#Backup
saveRDS(wrapped_gridIndoPacific, "Completeness_data/wrapped_gridIndoPacific.rds")
saveRDS(IndoPacific_SS, "Completeness_data/IndoPacific_SS_140324.rds")
saveRDS(IndoPacific_Cell_Assembly, "Completeness_data/IndoPacific_Cell_Assembly.rds")
#End of this codeEstimate sample-coverage completeness per sampling unit
Then, we estimated the completeness by sampling units (Cell id in a grid of ~\(100~ km^2\)). This was also conducted in the HiperGator (Do not run):
### Assess completeness of the sampling units
#* Orlando Acevedo-Charry
#* Updated on March 22nd, 2024
#Packages
library(tidyverse)
library(iNEXT) #rarefaction
#Load the data from Code_Islands_eBird_filtering.R
Caribbean_SS <- readRDS("Completeness_data/Caribbean_SS_140324.rds")
IndoPacific_SS <- readRDS("Completeness_data/IndoPacific_SS_140324.rds")
Caribbean_Cell_Assembly <- readRDS("Completeness_data/Caribbean_Cell_Assembly.rds")
IndoPacific_Cell_Assembly <- readRDS("Completeness_data/IndoPacific_Cell_Assembly.rds")
wrapped_gridCaribbean <- readRDS("Completeness_data/wrapped_gridCaribbean.rds")
wrapped_gridIndoPacific <- readRDS("Completeness_data/wrapped_gridIndoPacific.rds")
# Completeness with iNEXT and new grid ####
#Caribbean region
names(Caribbean_Cell_Assembly) #each row is a species (2254), and columns 2:10725 is the cell ID
listCaribbean<-as.data.frame(Caribbean_Cell_Assembly[,2:ncol(Caribbean_Cell_Assembly)])
#iNterpolation and EXTrapolation for the first three Hill numbers (see Chao et al. 2004):
Caribbean_iNEXT<-iNEXT(listCaribbean,q=c(0,1,2),
datatype = "abundance",
endpoint = NULL, #double the reference sample size
conf=0.95, nboot=500)
saveRDS(Caribbean_Completeness, "Completeness_output/Caribbean_Completeness.rds") #DONE
CompleteCaribbeanCells <- Caribbean_iNEXT$iNextEst$coverage_based %>%
filter(Method == "Observed" & Order.q %in% 0 & SC >= 0.9)
#Exclude cells in SS with less than 4 species
Caribbean_SS = Caribbean_SS[Caribbean_SS$cell %in% CompleteCaribbeanCells$Assemblage, ]
saveRDS(Caribbean_SS, "Completeness_output/Caribbean_SS_Completeness.rds") #DONE
#Get the grid cell boundaries for cells with completeness >0.9
wrapped_gridCaribbean = wrapped_gridCaribbean[wrapped_gridCaribbean$seqnum %in% CompleteCaribbeanCells$Assemblage, ]
gridCaribbean = merge(wrapped_gridCaribbean,CompleteCaribbeanCells,by.x="seqnum",by.y="Assemblage")
saveRDS(gridCaribbean, "Completeness_output/gridCaribbean.rds") #DONE
#Now with the IndoPacific
listIndoPacific<-as.data.frame(IndoPacific_Cell_Assembly[,2:ncol(IndoPacific_Cell_Assembly)])
#iNterpolation and EXTrapolation for the first three Hill numbers (see Chao et al. 2004):
IndoPacific_iNEXT<-iNEXT(listIndoPacific,q=c(0,1,2),
datatype = "abundance",
endpoint = NULL, #double the reference sample size
conf=0.95, nboot=500)
saveRDS(IndoPacific_iNEXT, "Completeness_output/IndoPacific_Completeness.rds") #DONE
CompleteIndoPacificCells <- IndoPacific_iNEXT$iNextEst$coverage_based %>%
filter(Method == "Observed" & Order.q %in% 0 & SC >= 0.9)
#Exclude cells in SS with less than 4 species
IndoPacific_SS = IndoPacific_SS[IndoPacific_SS$cell %in% CompleteIndoPacificCells$Assemblage, ]
saveRDS(IndoPacific_SS, "Completeness_output/IndoPacific_SS_Completeness.rds") #DONE
#Get the grid cell boundaries for cells with completeness >0.9
wrapped_gridIndoPacific = wrapped_gridIndoPacific[wrapped_gridIndoPacific$seqnum %in% CompleteIndoPacificCells$Assemblage, ]
gridIndoPacific = merge(wrapped_gridIndoPacific,CompleteIndoPacificCells,by.x="seqnum",by.y="Assemblage")
saveRDS(gridIndoPacific, "Completeness_output/gridIndoPacific.rds") #DONE
#End of this codeExample: Bananaquit populations
We selected a widespread distributed species to test our predictions of difference in stationary dynamic abundance from the mainland to islands, the Bananaquit Coereba flaveola
library(tidyverse) #seven packages in one for data management and visualization
library(MASS); #R functions and datasets to support "Modern Applied Statistics with S", a book from W.N. Venables and B.D. Ripley
source("ROUSSE-2.0.R") #Code from Dennis & Ponciano Ecology 2014.
#Call data and manipulate it
Coefla <- readRDS("Completeness_output/Caribbean_SS_Completeness_elevation.rds") |>
ungroup() |>
filter(scientific_name == "Coereba flaveola") |>
mutate(Day.t = as.numeric(observation_date)-15339, #First estimation of estimate in a year basis
Site.k = as.character(cell),
Year.t = as.numeric(year))
#Filter cells with more than 3 years of sampling
CellsMore3yrs <- Coefla %>%
group_by(Site.k, Year.t) %>%
summarise(Lists = n()) %>%
group_by(Site.k) %>%
summarise(List_Year = n()) %>%
filter(List_Year >= 3)
#Select only cells with more than 3 years of data
Coefla = Coefla[Coefla$Site.k %in% CellsMore3yrs$Site.k, ]
#Mean count per year per site
Coefla <- Coefla |>
group_by(Site.k, Year.t) |>
summarise(Observed.year.y = round(mean(observation_count),0)) |>
left_join(Coefla, by = c("Site.k","Year.t"))
#Mean count per day per site
Coefla <- Coefla |>
group_by(Site.k, Day.t) |>
summarise(Observed.day.y = round(mean(observation_count),0)) |>
left_join(Coefla, by = c("Site.k","Day.t"))
Coefla |>
dplyr::select(state_code, Site.k, Year.t, Day.t, Observed.year.y, Observed.day.y, observation_count) |>
arrange(state_code) |>
head()
saveRDS(Coefla, file = "Coereba_flaveola_data.rds")We can use the map of coast line Mainlan-Island map from Sayer et al. (J. Oper. Oceanogr.)
library(sf); #deal with shape files## Linking to GEOS 3.11.0, GDAL 3.5.3, PROJ 9.1.0; sf_use_s2() is TRUElibrary(tidyverse)## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.5.0 ✔ tibble 3.2.1
## ✔ lubridate 1.9.2 ✔ tidyr 1.3.0
## ✔ purrr 1.0.2## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errorslibrary(maps); #load maps##
## Attaching package: 'maps'
##
## The following object is masked from 'package:purrr':
##
## maplibrary(ggmagnify); # Zoom detail (inset figures)
Coefla <- readRDS("Coereba_flaveola_data.rds")
#The grid of cells
gridCaribbeanPost <- readRDS("Completeness_output/gridCaribbean.rds") #3823 Sites
#Only have the values in Coefla
gridCaribbeanPost = gridCaribbeanPost[gridCaribbeanPost$seqnum %in% Coefla$Site.k, ] #1033 Sites (cells k)
#A world map as reference
world1 <- sf::st_as_sf(map(database = 'world', plot = FALSE, fill = TRUE))
world1## Simple feature collection with 253 features and 1 field
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -180 ymin: -85.19218 xmax: 190.2708 ymax: 83.59961
## Geodetic CRS: +proj=longlat +ellps=clrk66 +no_defs +type=crs
## First 10 features:
## ID geom
## Aruba Aruba MULTIPOLYGON (((-69.89912 1...
## Afghanistan Afghanistan MULTIPOLYGON (((74.89131 37...
## Angola Angola MULTIPOLYGON (((23.9665 -10...
## Anguilla Anguilla MULTIPOLYGON (((-63.00122 1...
## Albania Albania MULTIPOLYGON (((20.06396 42...
## Finland Finland MULTIPOLYGON (((20.61133 60...
## Andorra Andorra MULTIPOLYGON (((1.706055 42...
## United Arab Emirates United Arab Emirates MULTIPOLYGON (((53.92783 24...
## Argentina Argentina MULTIPOLYGON (((-64.54916 -...
## Armenia Armenia MULTIPOLYGON (((45.55235 40...#Global Shoreline Vector
#coast <- st_read("USGSEsriWCMC_GlobalIslands_v3/v10/globalislandsfix.gdb")
#Different projections!
#Extract CRS of world1 map
world1_crs <- st_crs(world1)
#Transform the geodatabase to match the CRS projection of world1
#coastal <- st_transform(coast, world1_crs)
ggplot() +
geom_sf(data = world1)+
geom_sf(data=gridCaribbeanPost,
aes(color = SpRichness,
fill = SpRichness)) +
scale_color_gradient2(low = alpha("#4575b4", 0.95),
mid = alpha("#ffffbf", 0.95), midpoint = 73,
high = alpha("#d73027", 0.95))+
scale_fill_gradient2(low = alpha("#4575b4", 0.75),
mid = alpha("#ffffbf", 0.75), midpoint = 73,
high = alpha("#d73027", 0.75))+
geom_point(data = Coefla, aes(x = longitude, y = latitude), alpha = 0.5, size = 0.001)+
coord_sf(xlim = c(-91, -59.4),
ylim = c(8, 27)) +
scale_x_continuous(breaks = seq(from = -90, to = -60, by = 15))+
labs(y = "Latitude",
x = "Longitude",
tag = "",
title = "Neotropical region",
subtitle = "Caribbean (high eBird coverage-based completeness)",
color = "Species richness",
fill = "Species richness") +
theme_classic()+
theme(legend.position = "bottom",
legend.direction = "horizontal")+
guides(fill = F)+
geom_magnify(data = Coefla, aes(x = longitude, y = latitude),
from = c(xmin = -67.3, xmax = -65.6, ymin = 17.9, ymax = 18.6),
to = c(xmin = -70, xmax = -59, ymin = 22, ymax = 26.5),
shadow = TRUE) 
So, how is the time-series of the observed process?
gridCaribbeanPost$Site.k <- as.character(gridCaribbeanPost$seqnum)
Coefla <- Coefla |>
left_join(gridCaribbeanPost)## Joining with `by = join_by(Site.k, seqnum, SpRichness, SC, m, Method, Order.q,
## qD, qD.LCL, qD.UCL, geometry)`ggplot(Coefla, aes(x = observation_date, group = Site.k))+
geom_line(aes(y = log(Observed.day.y), color = SpRichness), alpha = 0.25, linetype = "dotted")+
geom_line(aes(y = log(Observed.year.y), color = SpRichness), alpha = 0.25, linetype = "dashed")+
scale_color_gradient2(low = alpha("#4575b4", 0.95),
mid = alpha("#ffffbf", 0.95), midpoint = 73,
high = alpha("#d73027", 0.95))+
geom_point(aes(y = log(Observed.day.y)), alpha = 0.01, size = 0.05, shape = 1)+
geom_point(aes(y = log(Observed.year.y)), alpha = 0.01, size = 0.1, shape = 2)+
scale_y_continuous(expand = c(0,0))+
scale_x_date(expand = c(0,0))+
labs(x = "Observation date",
y = "Log(mean observed)",
color = "Species richness")+
theme_classic()+
theme(legend.position = "bottom")
This is very messy an difficult to see, lets reduce to values around the Min., 1st Qu., Median, Mean, 3rd Qu., and Max. values of the Species richness gradient.
#Use some examples that represent sp richness: different species richness
summary(Coefla$SpRichness)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 5.0 37.0 53.0 119.5 225.0 383.0CoeflaSimp <- Coefla |>
filter(SpRichness %in% c(12,14,16:29, 31:68,72,
104,106,112,115,138,
155,156,160,174,179,186:188,196,
200:202,205:209,226:231,235,237,243:247,
250:253,261,270,278:284,293:299,
303:322,330:380)) |> #Select cells with variable richness and more than 200 observations
group_by(SpRichness) |>
sample_n(200) |>
as.data.frame()
ggplot(CoeflaSimp, aes(x = observation_date, group = Site.k))+
geom_line(aes(y = log(Observed.day.y), color = SpRichness), linetype = "dotted")+
geom_line(aes(y = log(Observed.year.y), color = SpRichness), linetype = "dashed")+
scale_color_gradient2(low = alpha("#4575b4", 0.95),
mid = alpha("#ffffbf", 0.95), midpoint = 73,
high = alpha("#d73027", 0.95))+
geom_point(aes(y = log(Observed.day.y)), size = 0.05, shape = 1)+
geom_point(aes(y = log(Observed.year.y)), size = 0.1, shape = 2)+
scale_y_continuous(expand = c(0,0))+
scale_x_date(expand = c(0,0))+
labs(x = "Observation date",
y = "Log(mean observed)",
color = "Species richness")+
theme_classic()+
theme(legend.position = "bottom")
In the previous figure, we assumed that the mean \(\mu_{t,k}\) of counts each year (or day in the future) \(t\) in each site \(k\) represents the stationary distribution of the stochastic density dependence in the population.
Look at the lower values of observed abundance in high richness sites, and higher estimates in low richness sites. Could we try to plug-in these time-series in the Gompertz state-state model? The problem is that this data is not uniform… so we need the diffusion process of Ornstein-Uhlenbeck process (Dennis & Ponciano, Ecology 2018).
Lele et al. (Ecol. Lett. 2007) provide the example of the GSS model as a stochastic, density-dependent model (based on Dennis et al., Ecol. Monogr. 2006). In this model, the latent variable of \(N_t\) as the true abundance at time \(t\) in the form: \[ N_t = N_{t-1}* {\rm exp} \left(a + b * {\rm ln}N_{t-1} + E_t \right) \] where \(a\) is the intrinsic growth rate (a constant) and \(b\) is a density dependence parameter (another constant), while the environmental stochasticity is \(E_t \sim {\rm Normal}(0,\sigma^2)\). The variance \(\sigma^2\) is the process noise or environmental variability of the system (Lele et al., Ecol. Lett. 2007), where higher values could serve to identify catastrophic stochasticity (Dennis et al., Ecol. Monogr. 1991).
Thus, we can define \(X_t = {\rm ln}(N_t)\), or the logarithm of the unobserved population abundances at time \(t\), \(t=(0,1,...,q)\). Applying the logarithm scale:
\[ X_t = X_{t-1} + a + (b*X_{t-1}) + E_t \] factorizing out, \[ X_t = a + X_{t-1}(1+b) + E_t \] and simplify the terms, \(c=1+b\) defines the density dependence parameter:
\[ X_t = a + cX_{t-1} + E_t \]
In the GSS model, the probability distribution of \(X_t\) is normal with mean and variance that change in function of time. We begin with a first value \(X_0\), that is plug-in in \(X_1 = a + c(X_0) +E_1\), then this \(X_1\) is plug-in in \(X_2 = a + c(X_1) + E_2 = a + c(a + c(X_0) +E_1) + E_2\), and so on… eventually, when \(X_t \rightarrow{\infty}\), \(X_t\) approaches a time-independent stationary distribution (carrying capacity):
\[ X_{\infty} \sim {\rm Normal} \biggl(\mu=\frac{a}{(1-c)}, {\rm var} = \frac{\sigma^2}{(1-c)^2} \biggr) \] Thus, we can assume that our \(X_0\) is part of that stationary distribution for the long-run time-series.
The observations \(Y_t\) are related to the true population abundances \(X_t\) by: \(Y_t = X_t + F_t\), where \(F_t \sim {\rm Normal}(0,\tau^2)\). Here, the variance \(\tau^2\) quantifies the amount of measurement error or observation error. We have to compute the MLE for \(\Theta = [a, c, \sigma, \tau]\). However, we can also plug the observations \(Y_t\) as a Poisson process (\(Y_t \sim {\rm Poisson}(N_t)\), or in log-scale \(Y_t \sim {\rm Poisson}({\rm ln}(X_t))\)), reducing a parameter estimation (\(\tau\)). However, the nature of the data from eBird require account for the variaton of observation error in the state-space model. In addition, this data is unequal in time intervals, requiring including the diffusion Ornstein-Uhlenbeck process (Dennis & Ponciano, Ecology 2018).
head(CoeflaSimp)## Site.k Day.t Observed.day.y Year.t Observed.year.y checklist_id common_name
## 1 743050 1961 3 2017 3 S37017271 Bananaquit
## 2 743050 3667 1 2022 7 G7739403 Bananaquit
## 3 743050 1962 2 2017 3 S37017210 Bananaquit
## 4 598530 1222 2 2015 2 S23327166 Bananaquit
## 5 602177 2913 2 2019 2 G4829124 Bananaquit
## 6 743780 3017 2 2020 5 S66777258 Bananaquit
## scientific_name observation_count country state_code locality_id latitude
## 1 Coereba flaveola 3 Guadeloupe GP- L5863751 16.17101
## 2 Coereba flaveola 1 Guadeloupe GP- L2357226 16.16975
## 3 Coereba flaveola 2 Guadeloupe GP- L5863741 16.17204
## 4 Coereba flaveola 2 Bahamas BS-EX L3625605 24.71089
## 5 Coereba flaveola 2 Bahamas BS-EX L831459 24.38150
## 6 Coereba flaveola 2 Guadeloupe GP- L2391989 16.17064
## longitude observation_date time_observations_started
## 1 -61.12085 2017-05-14 12.350000
## 2 -61.11943 2022-01-14 14.216667
## 3 -61.11568 2017-05-15 8.516667
## 4 -76.82254 2015-05-06 17.000000
## 5 -76.62037 2019-12-22 11.000000
## 6 -61.10903 2020-04-04 12.033333
## observer_id sampling_event_identifier
## 1 obsr307806 S37017271
## 2 obsr206146,obsr1979263,obsr2901391 S100736979,S101011323,S101137954
## 3 obsr307806 S37017210
## 4 obsr196942 S23327166
## 5 obsr1378293,obsr1474099 S62941819,S62990988
## 6 obsr307806 S66777258
## protocol_type duration_minutes effort_distance_km number_observers
## 1 Stationary 5 0.000 1
## 2 Traveling 72 3.800 3
## 3 Stationary 11 0.000 1
## 4 Traveling 30 0.322 1
## 5 Traveling 60 0.600 2
## 6 Stationary 5 0.000 1
## all_species_reported group_identifier hour_sampling year month week
## 1 TRUE <NA> 12 2017 5 20
## 2 TRUE G7739403 14 2022 1 2
## 3 TRUE <NA> 9 2017 5 20
## 4 TRUE <NA> 17 2015 5 18
## 5 TRUE G4829124 11 2019 12 51
## 6 TRUE <NA> 12 2020 4 14
## day_of_year cell mu_count ID Elevation_STR_90m seqnum SpRichness
## 1 134 743050 3.071429 1328296 3 743050 12
## 2 14 743050 3.071429 8632992 3 743050 12
## 3 135 743050 3.071429 1433352 4 743050 12
## 4 126 598530 3.833333 539992 6 598530 12
## 5 356 602177 3.457143 10708100 5 602177 12
## 6 95 743780 5.350000 4546662 4 743780 12
## SC m Method Order.q qD qD.LCL qD.UCL
## 1 1.0000000 34 Observed 0 12 4.174038 19.82596
## 2 1.0000000 34 Observed 0 12 4.174038 19.82596
## 3 1.0000000 34 Observed 0 12 4.174038 19.82596
## 4 0.9227250 25 Observed 0 12 2.481840 21.51816
## 5 1.0000000 33 Observed 0 12 5.755010 18.24499
## 6 0.9266436 26 Observed 0 12 2.673333 21.32667
## geometry
## 1 POLYGON ((-61.09389 16.2064...
## 2 POLYGON ((-61.09389 16.2064...
## 3 POLYGON ((-61.09389 16.2064...
## 4 POLYGON ((-76.7697 24.74715...
## 5 POLYGON ((-76.57627 24.4039...
## 6 POLYGON ((-61.01139 16.1688...CoeflaSimpAnual <- CoeflaSimp |>
dplyr::select(Site.k, Year.t, Observed.year.y, SpRichness, country) |>
unique() |>
filter(Site.k != "2737266")
table(CoeflaSimpAnual$Site.k, CoeflaSimpAnual$Year.t) #this shows that some cells have only 1 or less observations##
## 2012 2013 2014 2015 2016 2017 2018 2019 2020 2021 2022 2023
## 2704578 0 0 0 0 0 1 0 1 0 0 0 0
## 2716904 0 0 0 0 0 0 0 1 0 0 0 0
## 2718358 0 0 0 0 1 1 1 1 1 1 1 1
## 2721275 1 1 1 1 1 1 1 1 1 1 1 1
## 2722003 0 0 0 1 0 0 1 0 1 0 1 1
## 2722005 1 0 0 0 1 1 1 0 1 1 1 1
## 2722726 0 0 0 0 1 0 0 0 0 0 0 0
## 2722732 0 0 1 0 1 1 1 1 1 1 1 1
## 2722733 0 1 1 1 1 1 1 1 1 1 1 1
## 2722734 0 0 1 0 1 1 1 1 1 1 1 1
## 2723462 1 0 0 0 0 1 1 0 0 1 1 1
## 2723463 0 0 0 0 0 0 1 0 0 0 1 1
## 2736545 0 0 0 1 1 0 1 1 0 0 0 0
## 2737996 0 0 0 1 1 2 2 2 2 2 2 1
## 2738030 0 0 0 0 0 0 0 1 0 0 0 0
## 2738759 0 0 0 0 0 0 0 1 0 1 1 0
## 2740215 1 1 1 1 1 1 1 1 1 1 1 1
## 2740216 0 0 0 0 0 0 0 0 0 0 0 1
## 2741648 0 0 0 1 0 1 1 1 1 1 1 0
## 2742392 0 0 0 0 0 1 0 0 0 1 0 0
## 2743848 1 1 0 0 1 1 1 1 1 1 1 1
## 2746771 1 1 1 1 1 1 1 1 0 1 1 1
## 2746772 0 0 0 0 1 1 1 0 1 0 0 0
## 2753407 0 0 1 1 1 1 1 1 1 1 0 0
## 2754136 0 1 1 1 1 1 1 1 1 1 1 1
## 2754137 1 0 0 0 1 0 1 0 0 0 0 1
## 2754865 0 0 0 0 0 0 0 0 0 0 1 1
## 2754866 1 1 1 1 1 1 1 1 1 1 1 1
## 2754867 0 0 1 1 1 1 1 1 1 1 1 1
## 2754868 0 0 1 1 1 0 0 1 1 1 1 1
## 2755597 0 1 0 0 0 0 0 0 1 0 0 1
## 2756337 0 0 0 1 0 1 1 1 1 1 1 1
## 2756338 1 0 0 0 0 1 0 1 0 1 1 0
## 2757068 0 1 0 0 0 1 0 1 1 1 1 1
## 2757069 0 1 1 1 1 1 1 1 1 1 1 1
## 2757070 0 0 0 0 0 1 1 1 1 1 1 1
## 2803707 0 0 0 1 0 1 1 1 1 1 1 1
## 2808766 0 0 0 0 0 0 0 1 1 1 1 0
## 2808806 0 0 1 1 1 1 1 1 1 1 1 1
## 2809535 0 1 0 0 1 1 1 0 1 1 1 1
## 2810226 1 1 1 1 1 1 1 1 1 1 1 1
## 2810230 0 0 1 1 1 1 1 1 1 1 1 1
## 2810956 0 1 1 1 1 1 1 1 1 1 1 1
## 2813875 0 0 1 1 1 1 1 1 1 1 1 1
## 2813881 1 1 1 0 0 1 1 1 1 1 1 1
## 2814603 1 1 1 1 1 1 1 1 1 1 1 1
## 2814605 0 1 1 1 1 1 1 1 1 1 1 1
## 2816069 1 1 1 1 1 1 1 1 1 1 1 1
## 2816798 1 1 1 1 1 1 1 1 1 1 1 1
## 2816799 0 0 1 1 1 1 1 1 1 1 1 1
## 2817526 1 1 1 1 1 1 1 1 1 1 1 1
## 2818258 0 1 0 1 1 1 1 1 1 1 1 1
## 2818987 0 0 0 1 1 1 1 1 1 1 1 1
## 2819718 0 0 0 1 1 1 1 1 1 1 1 1
## 2820437 1 0 0 0 0 0 0 0 0 0 0 1
## 2820448 0 0 0 0 0 1 1 1 1 1 1 1
## 2821166 0 0 0 1 0 1 0 1 0 1 1 1
## 2821168 0 0 0 0 1 1 1 1 1 1 1 1
## 2822633 0 0 0 0 0 1 1 1 1 1 1 1
## 2823363 0 1 0 1 1 1 1 1 1 1 1 1
## 2824088 1 1 1 1 1 1 1 1 1 1 1 1
## 2824094 1 1 1 1 1 1 1 1 1 1 1 1
## 2825550 0 0 0 0 0 0 0 1 1 0 1 1
## 2826280 0 0 1 1 1 1 1 1 1 1 1 1
## 2827736 1 1 1 1 1 1 1 1 1 1 1 1
## 2827747 0 1 1 1 1 1 1 1 1 1 1 1
## 2828477 0 0 1 1 0 1 1 1 1 1 1 1
## 2829207 1 0 1 1 1 1 1 1 1 1 1 1
## 2829937 0 0 0 0 0 0 0 0 2 2 2 2
## 2832840 0 0 0 0 1 1 1 1 1 1 1 1
## 2832841 0 1 0 1 1 1 1 1 1 1 1 1
## 2833572 1 1 0 0 1 1 1 1 1 1 1 1
## 2833573 0 0 0 1 1 1 1 1 1 1 1 1
## 2834299 0 0 0 1 1 1 1 1 1 1 1 1
## 2834300 0 0 0 1 1 1 1 1 1 1 1 1
## 2834304 1 1 0 1 1 1 1 1 1 1 1 1
## 2834306 1 1 1 1 1 1 1 1 1 1 1 1
## 2834307 1 1 0 1 1 1 1 1 1 1 1 1
## 2835030 1 1 1 1 1 1 1 1 1 1 1 1
## 2835031 1 0 1 1 1 1 1 1 1 1 1 1
## 2835036 1 0 0 1 1 1 1 1 1 1 1 1
## 2835764 0 0 1 1 1 1 1 1 1 1 1 1
## 2836492 0 1 0 0 0 0 0 1 1 1 1 1
## 2836493 0 0 1 1 1 1 1 1 1 1 1 1
## 2837953 0 0 0 0 0 0 1 1 1 1 1 1
## 2840143 0 0 0 1 0 1 1 1 1 1 1 1
## 2843066 0 0 0 0 0 0 0 0 0 0 0 1
## 2843088 1 1 1 1 1 1 1 1 1 1 1 1
## 2845278 0 0 0 1 0 0 0 0 1 0 0 1
## 2846008 1 0 1 1 1 1 0 1 0 1 1 1
## 2847460 1 1 1 1 1 1 1 1 1 1 1 1
## 2847461 1 1 1 1 1 1 1 1 1 1 1 1
## 576642 0 1 0 1 1 0 0 0 0 0 0 0
## 577371 0 0 0 0 0 0 1 0 0 0 0 0
## 577372 0 0 1 0 1 1 1 1 0 0 1 0
## 577374 0 0 0 1 0 0 0 1 0 0 0 0
## 579571 0 0 1 1 0 1 1 0 0 0 0 1
## 580300 0 0 1 1 1 1 1 1 0 0 0 0
## 581030 0 0 0 0 1 0 1 0 0 0 0 0
## 581761 0 0 0 1 1 1 1 1 0 0 0 0
## 582490 0 1 0 0 1 1 1 1 0 0 0 1
## 583218 1 0 0 1 1 1 1 1 0 0 0 0
## 583946 0 0 0 0 0 0 0 0 1 0 0 0
## 583947 0 0 0 0 0 0 0 0 0 1 1 0
## 583948 0 0 0 0 0 1 0 0 1 0 1 0
## 584672 0 0 0 0 0 1 0 1 0 0 1 1
## 584675 0 0 0 0 0 0 0 1 0 0 0 0
## 584676 0 0 1 0 1 0 1 1 0 1 0 1
## 585397 0 0 0 0 0 1 1 1 0 1 1 1
## 585402 1 0 0 1 1 0 1 1 0 0 1 1
## 585404 0 0 1 0 0 0 0 0 1 0 1 0
## 586133 0 0 0 0 1 0 0 0 0 0 0 0
## 590499 0 0 0 0 0 1 0 0 0 0 0 0
## 591228 0 1 0 0 1 1 0 0 0 0 0 0
## 591970 0 0 0 0 1 0 0 0 0 0 0 0
## 592686 0 1 0 0 0 0 1 0 0 0 0 0
## 592700 0 0 0 0 0 0 0 0 0 1 1 1
## 593416 0 0 0 0 0 0 0 0 0 1 0 1
## 593420 0 0 0 0 0 0 0 0 1 0 0 0
## 593421 0 0 0 0 0 0 0 1 0 0 0 0
## 594149 0 0 1 0 1 0 0 1 1 1 0 0
## 594150 0 0 1 0 1 0 0 0 0 1 1 1
## 594151 0 1 1 0 1 1 1 1 1 1 1 1
## 594152 0 0 0 0 0 0 0 1 0 1 0 0
## 594159 0 0 0 1 0 0 0 1 0 0 0 1
## 594875 0 0 0 0 0 0 0 0 0 0 1 1
## 594880 0 0 0 0 0 0 0 0 1 0 0 0
## 594881 0 0 0 0 0 0 0 1 0 0 0 0
## 594889 0 1 0 0 0 0 1 0 0 1 1 1
## 595604 0 0 0 0 0 0 0 0 0 0 1 0
## 595605 0 0 0 0 0 0 0 0 1 0 0 0
## 596349 0 0 0 0 0 1 0 0 0 0 0 0
## 597078 0 0 0 1 0 1 0 1 1 0 0 0
## 597079 0 1 0 0 0 0 1 1 1 1 1 0
## 598530 0 0 0 1 0 0 1 0 0 0 0 1
## 598537 0 0 0 1 0 1 0 0 0 0 1 1
## 599265 0 0 0 0 1 0 1 1 1 1 1 1
## 599266 0 0 0 1 0 0 0 0 1 0 0 0
## 599980 0 0 0 0 0 0 1 0 0 0 1 0
## 600709 0 0 0 0 1 0 0 0 0 0 0 0
## 601439 0 0 0 0 0 0 0 0 0 1 1 1
## 601452 0 1 0 1 0 0 1 1 0 0 1 0
## 602177 1 0 0 1 1 0 1 1 1 1 1 1
## 602182 0 0 0 0 0 1 0 0 0 0 0 1
## 602913 0 1 1 1 0 1 1 1 0 0 1 1
## 602915 0 0 0 0 0 0 0 1 0 0 0 1
## 605095 1 0 0 0 0 0 0 0 0 1 1 1
## 610932 0 0 0 0 0 0 1 1 0 0 1 1
## 611675 0 1 0 1 1 1 1 1 1 0 1 1
## 611676 0 1 0 1 1 0 1 1 1 1 1 1
## 612404 0 1 0 0 0 0 1 1 1 0 1 0
## 613121 0 0 0 0 0 0 0 1 1 0 0 1
## 613133 1 0 0 1 0 0 1 1 1 1 1 1
## 613134 0 1 0 0 0 1 1 1 1 0 1 0
## 618961 0 0 1 1 0 0 1 1 0 0 1 0
## 638671 0 0 0 1 0 0 0 1 0 1 1 1
## 638674 0 0 0 1 0 0 0 0 0 0 1 0
## 639401 1 1 0 0 0 1 0 1 1 1 1 1
## 639404 0 1 1 1 1 1 0 1 0 1 1 1
## 640130 0 1 0 0 0 1 0 0 0 1 1 1
## 640131 0 0 1 1 1 1 1 1 0 1 1 1
## 640863 0 0 1 1 0 0 1 1 1 0 1 1
## 641574 0 0 0 1 0 1 0 0 0 0 0 0
## 641575 0 0 0 0 1 1 1 0 0 0 0 0
## 642303 0 0 0 1 1 1 1 1 1 1 1 1
## 642306 0 0 0 1 0 1 0 0 0 0 0 0
## 643033 0 0 0 1 1 1 0 1 0 1 1 0
## 647358 0 0 0 0 0 1 1 0 0 1 1 1
## 648087 0 1 1 1 0 1 1 1 1 0 1 1
## 648088 1 0 0 0 0 0 0 0 1 0 1 0
## 648090 0 0 0 0 0 0 1 0 0 1 0 0
## 648091 0 0 0 0 0 0 0 1 1 0 1 0
## 648092 1 0 0 0 1 1 1 1 1 1 1 1
## 648817 0 0 0 0 0 0 1 0 0 0 0 0
## 648820 0 1 0 0 0 0 0 0 0 0 0 0
## 648821 1 1 0 1 1 1 1 1 0 1 1 1
## 648823 0 0 0 0 0 0 0 0 0 0 1 0
## 649553 0 0 1 0 1 0 0 0 1 0 1 1
## 649554 0 0 0 0 0 1 1 1 1 0 0 0
## 650278 0 0 0 0 0 0 0 1 0 0 0 0
## 650285 0 0 1 1 0 0 1 1 0 1 1 1
## 651007 0 1 1 0 1 0 1 0 1 1 1 1
## 651010 0 0 0 1 0 1 0 0 0 0 0 0
## 651013 0 0 0 1 1 1 1 1 1 0 0 1
## 651014 0 0 0 0 0 0 0 0 1 1 1 0
## 651015 0 0 0 0 1 0 0 0 0 0 1 0
## 651737 0 1 0 0 0 1 0 0 1 0 0 0
## 651739 0 0 0 0 1 0 0 0 0 0 0 0
## 651745 0 0 0 0 0 1 1 1 0 0 0 0
## 651746 0 1 1 0 0 0 1 0 0 1 1 0
## 652467 0 0 0 0 0 1 0 0 1 0 0 0
## 652468 0 0 1 0 0 0 0 1 1 0 0 0
## 652475 0 0 0 0 0 0 0 1 1 0 0 1
## 653199 0 1 0 1 1 0 0 0 1 0 1 1
## 653936 0 0 0 1 1 1 0 1 1 1 1 0
## 654665 0 0 0 0 0 0 0 0 1 0 0 0
## 655394 0 0 0 0 0 1 0 1 0 0 1 1
## 655395 0 0 0 0 0 0 0 1 0 0 0 0
## 656120 0 0 0 0 0 0 0 0 1 0 0 1
## 656122 0 0 0 0 0 0 0 0 0 0 1 0
## 656123 1 1 0 1 1 1 1 1 1 1 1 1
## 656124 0 0 1 1 1 1 1 1 1 1 1 1
## 656126 0 0 0 1 0 0 0 0 0 1 1 1
## 656850 0 0 0 0 0 0 0 1 0 0 0 0
## 656851 0 0 0 0 0 0 1 0 0 0 1 0
## 656852 0 0 0 1 0 1 1 0 1 1 1 1
## 656853 0 0 0 0 0 0 0 1 1 1 1 1
## 656854 0 0 0 0 0 0 0 0 1 0 0 0
## 656856 0 1 1 0 0 1 0 1 0 0 1 1
## 657577 0 0 1 0 0 1 0 1 0 0 0 0
## 657580 0 0 0 0 0 0 0 0 0 0 0 1
## 657585 0 0 0 1 0 1 1 1 1 0 0 1
## 657586 0 0 0 0 0 0 1 0 0 0 1 1
## 657626 0 0 0 0 0 0 1 1 1 0 1 0
## 657627 0 0 0 0 0 1 0 1 0 0 1 1
## 657628 0 0 0 0 0 0 0 1 0 0 0 0
## 658355 0 0 0 0 0 0 1 0 0 0 0 0
## 658356 0 0 0 0 0 0 1 1 1 0 1 0
## 659043 0 0 0 0 0 0 0 0 0 0 0 1
## 659086 0 0 0 0 0 0 0 1 0 0 0 0
## 659091 0 0 0 0 0 0 1 0 0 1 0 1
## 659820 0 0 0 1 0 0 0 0 1 0 0 0
## 659821 0 0 0 0 0 1 0 0 0 0 0 0
## 660549 0 0 1 1 0 0 1 0 0 0 0 1
## 660554 0 0 1 0 0 0 0 0 0 0 0 0
## 661285 0 0 0 0 0 1 0 0 1 0 0 0
## 662017 0 0 0 0 0 0 1 0 0 0 1 0
## 662746 1 0 0 1 0 1 0 0 0 1 1 1
## 663455 0 0 0 0 0 0 0 0 0 1 0 0
## 663476 0 1 0 0 0 0 1 1 1 0 1 1
## 663477 0 0 1 0 1 1 1 1 1 1 0 1
## 663478 0 0 1 0 0 0 0 1 1 1 1 0
## 664207 0 0 0 0 0 0 0 0 0 1 1 0
## 664933 0 0 0 0 0 0 0 0 0 1 0 0
## 664936 0 0 0 0 0 0 1 0 0 1 0 1
## 665632 0 0 0 0 0 0 1 0 0 0 0 0
## 665662 1 1 0 0 1 0 0 0 0 1 0 1
## 665663 0 0 0 0 0 0 0 1 0 0 0 1
## 665668 0 0 0 0 0 0 0 0 0 1 0 0
## 666361 1 1 1 1 1 1 1 1 1 0 1 0
## 666392 0 0 0 0 0 0 0 0 1 1 0 1
## 667093 0 0 1 1 1 0 0 1 0 0 1 0
## 667123 0 0 0 0 0 0 0 0 0 1 0 1
## 667125 0 1 0 0 0 0 0 0 0 0 0 1
## 667128 0 0 0 0 0 0 0 0 0 1 1 0
## 667854 0 0 0 0 0 0 0 1 1 1 1 1
## 668563 1 0 1 1 1 1 1 1 1 1 1 0
## 668577 0 0 0 0 0 0 0 0 0 0 0 1
## 668579 0 0 0 0 0 0 0 0 0 0 0 1
## 668584 0 0 0 0 1 0 0 0 0 1 0 0
## 668587 0 0 0 0 0 0 0 0 0 0 1 0
## 669289 0 0 0 0 1 1 1 1 0 0 0 0
## 669292 0 1 0 1 0 1 1 1 1 0 0 0
## 669297 0 0 0 0 0 0 0 0 0 1 0 0
## 669308 0 0 0 0 0 1 1 1 0 0 0 0
## 669309 0 0 0 1 1 1 0 1 1 0 1 1
## 670021 1 0 0 0 1 0 1 1 1 1 0 0
## 670037 0 0 0 0 0 0 0 1 1 0 0 0
## 670038 0 0 0 0 0 0 0 0 0 1 0 0
## 670039 0 0 0 0 0 0 0 0 0 0 1 0
## 670040 0 0 0 0 0 0 0 0 0 0 0 1
## 670048 0 0 0 0 0 0 0 0 1 0 0 1
## 670748 0 0 0 0 1 1 1 0 0 0 0 0
## 670749 1 0 0 0 0 0 1 0 0 0 0 0
## 670750 0 0 1 0 1 1 1 1 0 1 1 0
## 670766 0 0 0 0 0 0 0 0 0 0 1 0
## 670776 0 0 0 0 0 0 1 0 0 0 0 0
## 670778 1 0 0 0 1 0 0 1 0 1 0 0
## 671495 0 0 0 1 0 0 0 0 1 1 0 0
## 671497 0 0 0 0 0 0 0 0 1 0 1 0
## 672214 0 0 0 0 0 0 0 0 0 1 0 0
## 672215 0 0 0 0 0 0 0 0 0 0 0 1
## 672227 0 0 0 0 0 0 1 0 1 1 0 0
## 672237 0 0 0 0 0 0 1 0 0 0 0 0
## 672238 0 0 0 0 0 0 0 0 1 0 0 0
## 672942 0 0 0 0 0 0 0 1 0 1 1 0
## 672943 1 1 1 0 0 1 0 0 0 1 1 1
## 672944 1 1 1 1 1 1 1 1 1 1 1 1
## 672945 1 0 0 0 0 0 0 0 0 0 1 0
## 672965 0 1 1 0 1 1 1 1 1 1 1 1
## 672966 0 0 0 0 0 0 0 0 0 1 0 1
## 673672 0 0 0 0 0 0 0 0 0 0 0 1
## 674401 0 0 0 0 0 0 0 1 0 0 0 0
## 674402 0 0 0 1 0 0 0 0 0 0 1 0
## 674405 0 1 0 1 0 0 0 0 0 0 1 0
## 674406 0 1 0 0 0 0 0 0 0 0 0 0
## 674420 0 0 1 0 0 0 0 0 0 0 1 1
## 675130 0 1 0 0 0 1 1 0 0 0 1 1
## 675131 0 0 0 0 0 0 0 0 1 0 0 0
## 675864 0 0 0 0 0 0 0 0 1 1 1 1
## 675877 0 0 0 0 0 0 0 1 1 1 0 0
## 676589 0 0 0 0 0 1 0 0 0 0 0 0
## 676593 0 0 0 0 0 0 0 0 0 0 1 0
## 676594 0 1 0 0 0 0 0 0 1 0 1 0
## 676605 0 0 0 0 0 0 0 0 0 0 0 1
## 676606 1 1 1 1 1 0 0 0 1 1 1 1
## 676608 0 0 0 0 0 1 0 0 0 0 0 0
## 677323 0 0 0 0 0 0 0 1 1 0 0 0
## 677329 0 0 0 0 0 0 0 0 1 1 1 1
## 677334 0 0 0 0 0 0 0 0 1 1 0 0
## 677335 0 0 0 0 0 0 1 1 1 1 1 1
## 677336 1 0 1 1 0 1 1 1 1 0 1 1
## 678049 0 0 0 0 0 0 0 0 1 0 0 0
## 678779 0 0 0 0 0 1 0 0 0 0 0 0
## 678801 0 0 0 1 0 0 0 0 0 0 0 0
## 679530 0 0 1 0 0 0 0 0 0 0 0 0
## 680261 0 0 0 1 0 0 0 0 1 0 0 1
## 680266 0 0 0 0 0 0 0 0 0 0 1 1
## 680991 0 0 0 1 1 0 0 0 1 0 1 1
## 680996 0 0 0 0 0 0 0 0 0 0 0 1
## 681720 0 1 0 0 0 1 1 1 1 1 1 1
## 681725 0 1 1 1 1 1 1 1 1 1 1 1
## 682455 0 1 0 1 1 1 1 1 1 1 1 1
## 683180 0 0 0 0 0 0 1 1 0 1 0 1
## 683181 1 0 0 0 0 0 0 0 0 0 0 0
## 683184 1 1 1 1 1 1 1 1 1 1 1 1
## 683909 0 0 0 0 0 1 0 0 0 0 0 0
## 691218 0 0 1 0 0 1 1 1 1 1 1 1
## 691946 0 1 1 1 1 1 1 1 1 1 1 1
## 691947 0 0 0 0 0 0 0 0 0 0 1 0
## 691948 1 0 1 0 1 1 1 1 1 1 1 1
## 691949 0 0 0 1 1 0 0 1 1 1 1 0
## 692676 1 1 1 1 1 1 1 1 1 0 1 0
## 692678 1 1 1 0 1 1 1 1 1 0 0 0
## 692679 0 0 0 0 1 1 1 1 1 1 1 1
## 693405 0 0 1 1 0 1 1 1 1 1 1 1
## 693406 0 1 0 0 0 0 1 1 1 1 1 0
## 693407 1 0 0 1 1 0 1 1 1 1 1 0
## 693408 0 1 0 0 0 1 1 1 1 1 1 1
## 693409 0 1 1 1 1 1 1 1 1 1 1 1
## 694133 0 1 0 0 1 0 1 0 1 1 1 0
## 694134 1 0 0 1 1 0 0 1 1 0 0 0
## 694135 0 0 0 0 0 0 0 1 0 1 1 0
## 694136 0 0 0 0 1 0 0 1 1 1 1 1
## 694137 0 0 0 0 0 0 0 1 1 0 0 1
## 694138 0 0 0 0 0 0 0 1 0 1 1 0
## 694139 1 0 0 0 0 1 1 1 1 1 1 1
## 694140 1 0 0 0 1 0 0 1 1 1 1 1
## 694862 0 0 1 1 1 0 0 1 1 1 1 1
## 694863 1 1 1 1 1 1 1 1 1 1 1 1
## 694864 0 0 0 0 0 0 0 0 0 1 0 1
## 694865 1 1 1 1 1 1 1 1 1 1 1 1
## 694866 0 0 0 1 0 0 1 1 0 1 1 1
## 694867 0 1 0 1 1 0 0 0 0 1 0 1
## 694868 0 0 0 0 0 0 0 0 0 1 0 0
## 694869 0 1 1 1 1 1 1 1 1 0 1 0
## 694870 0 1 1 1 1 1 1 0 1 1 1 1
## 695592 0 0 1 0 0 0 0 0 1 1 0 1
## 695593 1 1 1 1 1 1 1 1 1 1 1 1
## 695594 0 0 0 1 1 0 0 1 1 0 1 1
## 695595 0 1 0 0 1 0 1 0 0 1 0 0
## 695596 0 0 1 1 0 1 0 0 0 1 0 0
## 695597 0 0 0 1 0 0 0 0 0 0 0 0
## 695598 1 0 1 1 1 1 1 1 1 1 1 1
## 695599 0 0 0 0 1 0 1 1 1 1 1 0
## 695600 1 1 1 1 1 1 1 1 1 1 1 1
## 695601 0 1 1 0 0 0 0 1 1 0 1 1
## 696323 0 1 0 0 1 0 1 0 0 0 1 1
## 696325 1 0 0 0 1 0 0 0 1 0 0 0
## 696326 0 0 0 0 1 0 1 1 1 1 1 1
## 696327 0 0 0 0 0 0 0 1 1 0 1 1
## 696328 1 0 1 1 1 1 1 0 0 1 1 1
## 696329 1 1 1 1 0 0 1 1 1 1 1 0
## 696330 1 1 1 1 0 1 1 1 1 1 1 1
## 696331 0 1 0 0 0 1 0 1 0 1 1 1
## 697053 0 0 0 1 0 0 0 1 0 1 1 0
## 697054 1 1 1 1 1 1 1 1 1 1 1 1
## 697055 1 1 1 1 1 1 0 0 1 1 0 0
## 697056 1 1 1 1 1 0 1 0 0 1 1 1
## 697057 0 0 0 0 0 0 0 1 0 0 1 1
## 697058 0 0 0 1 0 0 0 1 0 0 0 0
## 697059 0 1 1 0 0 1 0 1 0 1 1 1
## 697060 0 0 0 0 0 0 0 0 1 1 1 0
## 697061 1 0 1 1 1 1 1 1 1 1 1 1
## 697784 0 1 1 0 0 0 0 0 0 0 0 0
## 697785 0 1 1 0 0 1 1 1 1 1 1 1
## 697786 0 1 1 0 1 1 0 0 1 1 1 1
## 697787 0 0 0 0 0 0 0 0 1 0 0 1
## 697788 0 0 0 1 0 0 0 1 1 0 1 0
## 697789 1 0 0 0 0 0 0 0 0 0 1 0
## 697790 0 0 0 0 1 0 0 0 0 0 1 0
## 697791 0 0 0 0 1 0 1 1 1 1 1 1
## 697792 1 1 1 1 1 1 1 1 1 1 1 1
## 698515 0 1 0 0 1 1 0 1 1 0 1 1
## 698516 0 0 0 0 1 1 1 0 1 1 0 1
## 698517 0 0 0 0 0 0 1 0 0 0 0 1
## 698518 1 0 0 1 1 0 1 1 0 1 1 0
## 698519 0 0 0 1 1 0 0 1 0 1 1 1
## 698520 0 0 0 0 0 1 0 0 1 1 1 1
## 698521 1 1 1 1 1 1 1 1 1 1 1 1
## 698522 0 1 0 0 0 1 1 1 1 1 1 1
## 699245 0 0 0 0 0 0 0 0 1 0 0 1
## 699247 0 0 1 0 0 0 0 1 1 1 1 0
## 699248 0 0 0 1 1 0 0 0 1 1 1 1
## 699249 1 0 1 1 1 0 1 1 0 1 1 1
## 699250 0 0 0 0 0 0 0 1 1 1 1 1
## 699251 1 1 1 1 0 1 1 1 1 1 1 1
## 699252 1 1 1 1 1 1 1 1 1 1 1 1
## 699976 1 0 0 0 0 0 1 0 0 0 0 0
## 699977 1 0 0 0 0 0 1 0 0 0 1 0
## 699978 0 0 0 0 0 1 0 0 1 1 1 1
## 699979 0 0 0 1 0 0 0 0 0 1 0 1
## 699980 0 1 0 0 1 1 1 1 1 1 1 1
## 699981 0 0 0 0 0 0 1 0 0 0 1 1
## 699982 1 1 0 0 1 0 1 1 1 1 1 1
## 699983 0 0 0 0 0 0 0 0 0 0 0 1
## 700706 0 0 0 0 1 0 1 1 1 1 1 1
## 700707 0 0 0 0 0 0 1 0 0 1 1 1
## 700708 0 1 1 0 0 1 1 1 1 1 1 1
## 700709 0 0 0 0 0 0 1 0 1 1 1 1
## 700710 0 0 0 1 0 0 0 1 0 1 1 1
## 700711 0 0 0 1 1 1 1 0 1 1 1 1
## 700712 1 1 1 1 1 1 1 1 1 1 1 1
## 700713 0 0 0 0 0 0 0 0 0 1 1 1
## 701436 1 1 1 1 1 1 1 1 1 1 1 1
## 701437 0 0 0 0 0 1 1 1 0 1 1 1
## 701438 0 0 0 0 1 0 0 0 1 1 1 1
## 701439 1 0 0 0 1 0 1 0 0 0 1 0
## 701440 0 0 0 0 0 0 0 1 1 0 1 0
## 701441 0 1 1 1 1 1 1 1 1 1 1 1
## 701442 0 0 0 1 0 1 1 1 0 1 1 1
## 701443 0 0 0 0 1 0 0 1 1 0 1 1
## 702167 0 0 0 0 0 0 0 0 0 1 0 0
## 702169 0 0 0 0 0 0 0 0 0 1 1 0
## 702170 1 1 1 1 1 1 1 1 1 1 1 1
## 702171 0 0 0 0 0 0 0 0 0 1 1 0
## 702172 0 1 1 0 1 1 0 1 1 1 1 1
## 702898 0 1 1 1 0 0 0 0 1 1 1 0
## 702899 0 0 0 1 0 0 0 0 0 0 0 1
## 702901 0 0 0 0 0 1 0 1 1 1 1 1
## 703630 1 0 0 0 0 1 0 0 1 1 1 1
## 703631 0 0 0 1 0 1 0 0 0 1 1 1
## 703633 0 1 0 0 0 0 1 0 0 0 0 0
## 703634 0 0 0 1 1 1 1 1 0 0 1 1
## 704360 0 0 0 0 0 0 0 0 0 1 1 1
## 704361 1 1 1 0 1 1 0 1 1 1 1 1
## 704363 0 0 1 1 0 1 1 1 0 0 1 1
## 705091 0 0 0 1 0 0 0 1 0 1 1 0
## 705095 0 0 0 0 0 0 0 1 0 0 0 0
## 705825 1 1 1 1 1 1 1 1 1 1 1 1
## 705826 0 0 0 0 0 0 0 1 1 1 1 1
## 705827 0 0 0 0 1 0 0 1 0 0 1 1
## 705832 0 0 0 1 0 1 0 0 0 0 1 1
## 706555 0 1 1 1 1 0 0 1 0 1 1 1
## 706556 1 1 1 1 1 1 1 1 1 1 1 1
## 706557 1 0 1 0 0 0 0 0 0 0 0 0
## 707286 1 1 1 1 1 1 1 1 1 1 1 1
## 707287 1 0 1 0 1 0 0 0 0 0 1 1
## 707288 1 1 0 0 1 1 0 0 0 0 1 1
## 708017 0 0 0 0 0 0 0 0 0 1 0 1
## 708018 0 1 0 1 1 0 0 1 1 0 1 1
## 708019 1 0 0 0 0 0 0 1 0 0 0 0
## 710200 0 0 0 0 0 0 1 1 0 0 0 1
## 710929 0 0 0 1 1 1 1 1 1 1 1 1
## 710930 0 0 0 1 1 1 1 1 1 1 1 1
## 711659 0 1 1 1 1 1 0 0 1 0 0 1
## 711660 1 1 1 0 1 1 1 1 1 1 1 1
## 711661 1 0 0 1 1 1 1 1 1 1 1 1
## 712390 0 0 1 0 0 1 1 1 1 1 1 1
## 717511 1 1 1 1 1 1 1 1 1 1 1 1
## 717512 0 1 1 1 1 1 1 1 1 1 1 1
## 718240 0 0 1 0 0 1 1 2 0 2 1 2
## 718241 0 1 1 1 1 1 0 1 1 1 1 1
## 718970 1 1 2 1 1 1 1 2 1 1 2 2
## 720425 0 0 0 0 1 0 0 1 0 0 0 0
## 720430 0 0 0 0 0 1 0 1 1 1 1 1
## 721155 0 1 0 0 1 0 0 0 1 0 0 0
## 721160 0 0 0 0 0 0 0 0 1 1 1 1
## 723345 0 0 1 0 1 1 1 0 1 1 1 0
## 723942 0 0 0 0 1 1 1 1 1 0 1 1
## 724679 0 0 1 1 0 1 1 1 1 1 1 1
## 724680 0 0 0 0 1 1 1 1 0 1 1 1
## 724703 0 0 0 0 0 0 0 0 1 0 0 0
## 725360 0 0 1 1 1 1 1 1 1 1 1 1
## 725534 0 0 0 0 0 1 1 1 0 1 1 1
## 726124 0 0 0 0 0 1 1 0 0 1 1 1
## 726137 0 0 0 0 1 0 0 0 0 0 0 1
## 726264 0 1 1 1 1 0 1 1 1 1 1 1
## 726853 0 0 0 0 0 0 0 1 0 1 1 1
## 726873 0 0 0 0 0 1 1 1 1 1 1 1
## 726878 1 1 1 1 1 1 1 1 1 1 1 1
## 726885 0 0 0 0 0 0 0 1 0 0 0 0
## 726994 0 1 1 0 1 0 0 1 0 1 1 1
## 727604 0 0 0 0 0 1 0 1 0 0 1 0
## 727605 0 1 0 0 0 0 0 0 0 0 0 1
## 727609 0 0 0 0 0 0 0 0 0 0 0 1
## 727626 0 0 0 0 0 0 1 1 1 0 0 0
## 727723 0 1 0 0 1 0 0 1 1 1 0 1
## 728453 0 0 0 0 1 1 1 1 0 1 0 0
## 728461 0 0 0 0 0 0 0 0 0 1 0 0
## 728462 0 1 0 0 0 1 0 0 1 0 1 1
## 729098 1 1 1 1 1 1 1 1 1 1 1 1
## 729191 0 0 0 0 0 0 1 0 1 0 1 1
## 729827 0 0 1 1 1 1 1 1 0 1 1 1
## 729828 0 1 0 0 0 1 1 1 1 1 1 1
## 730557 0 0 1 1 1 1 1 1 1 1 1 1
## 730558 0 0 0 0 0 0 0 0 1 0 1 1
## 731287 0 0 0 0 0 0 0 0 0 1 1 0
## 732106 0 0 1 0 1 0 0 1 1 1 1 1
## 732107 0 1 0 0 1 0 0 1 1 1 1 1
## 732831 0 0 1 1 1 1 1 0 1 1 1 1
## 732835 1 0 0 0 1 0 1 1 0 1 1 0
## 732836 0 0 0 0 0 1 1 0 1 1 1 1
## 732837 0 0 0 0 0 1 1 1 0 0 0 1
## 733560 0 0 0 1 1 1 0 0 0 0 0 0
## 733565 0 0 0 0 0 0 1 1 1 1 1 1
## 733566 0 0 1 1 0 0 1 1 1 1 1 1
## 733567 0 0 0 0 1 0 0 1 0 1 1 1
## 734211 0 0 0 1 1 1 1 1 1 1 1 1
## 734928 0 0 0 0 0 0 0 0 0 0 1 0
## 734931 0 0 0 0 0 1 0 0 0 1 1 1
## 734940 0 0 0 0 0 0 0 1 1 0 0 0
## 735590 1 0 0 0 0 0 0 0 0 0 0 0
## 735670 1 0 0 0 0 0 0 1 0 0 0 1
## 736387 0 0 0 0 0 0 0 0 0 0 1 0
## 736390 0 0 0 0 0 0 0 0 0 0 1 0
## 736399 0 0 0 0 0 0 1 1 1 1 1 1
## 736400 0 0 0 0 0 0 0 1 1 0 1 1
## 737129 1 0 0 1 1 1 1 1 1 1 1 1
## 737130 0 0 0 0 1 0 1 1 1 1 1 0
## 737782 0 0 0 0 0 0 0 0 1 0 1 0
## 737829 0 0 1 0 0 0 0 0 0 0 0 0
## 737859 0 0 0 1 1 1 1 1 1 1 1 1
## 737860 0 0 0 0 0 0 0 1 1 1 1 0
## 737864 0 1 1 0 1 1 1 1 0 1 1 1
## 738558 0 0 0 0 0 1 1 0 0 0 0 0
## 738589 0 0 0 0 0 0 0 1 0 1 1 0
## 738593 0 0 0 0 0 0 0 0 1 1 1 1
## 738594 0 1 1 1 0 1 1 0 1 1 1 1
## 738668 0 0 1 1 1 1 1 1 1 1 1 1
## 738669 0 0 0 0 1 1 0 1 1 1 1 1
## 738671 0 0 0 0 1 0 0 0 0 1 1 0
## 739289 0 0 0 0 0 0 0 0 1 1 1 1
## 739323 0 0 1 1 1 1 1 1 1 1 1 1
## 739324 0 0 0 0 0 0 0 0 1 1 1 1
## 739397 0 0 0 0 0 0 0 0 1 1 0 1
## 739398 0 0 1 1 0 1 1 0 1 0 0 0
## 739399 0 0 0 0 1 0 0 0 1 1 1 0
## 739400 0 0 0 1 1 1 0 0 0 1 0 0
## 739401 0 0 0 1 1 1 1 0 1 1 1 1
## 739402 1 0 0 0 0 0 0 0 0 0 1 1
## 740017 0 0 0 0 0 0 0 0 1 1 1 0
## 740052 0 0 0 0 1 0 1 1 0 0 1 1
## 740053 0 0 0 0 1 1 1 1 1 1 1 1
## 740127 0 1 0 0 1 1 1 1 1 1 0 1
## 740128 0 1 1 1 1 1 1 1 1 1 0 1
## 740129 0 0 1 0 1 0 0 0 1 0 0 1
## 740130 0 1 1 1 1 1 1 1 1 1 1 1
## 740131 0 0 0 0 1 0 0 0 0 1 0 0
## 740767 0 0 0 0 0 0 0 1 1 1 1 1
## 740768 0 0 0 0 0 0 0 0 1 1 1 1
## 740856 0 1 1 1 0 1 1 1 1 1 1 1
## 740857 0 0 0 0 0 0 0 0 0 1 1 1
## 740858 0 0 0 0 1 0 1 1 1 0 1 1
## 740859 0 1 1 1 1 0 1 1 1 1 1 1
## 740860 0 0 1 1 0 1 0 1 0 1 1 0
## 740861 0 0 0 1 0 0 0 1 0 0 1 0
## 741586 0 0 1 1 0 0 1 1 1 1 1 1
## 741587 0 0 0 0 0 0 0 1 1 0 1 0
## 741588 0 0 0 0 0 1 0 0 1 0 1 1
## 741590 0 0 0 1 1 1 0 0 0 0 1 1
## 742214 0 0 0 0 0 0 0 0 1 0 1 0
## 742227 0 0 0 0 0 0 0 1 1 1 1 0
## 742315 0 0 0 1 0 0 1 0 0 0 0 0
## 742320 0 0 1 1 1 0 1 1 1 0 1 1
## 742321 0 0 0 0 1 1 0 1 0 0 1 0
## 742322 0 0 1 1 1 1 1 0 1 1 0 0
## 743044 0 0 0 0 0 0 1 1 1 1 0 0
## 743045 0 0 0 0 1 0 1 0 0 0 0 1
## 743050 0 0 1 1 1 1 1 1 1 1 1 0
## 743051 0 1 1 1 0 0 0 0 1 1 1 0
## 743052 0 0 0 1 1 1 1 1 1 1 0 0
## 743777 0 0 0 0 0 0 1 1 1 1 1 1
## 743780 0 1 1 1 1 1 1 1 1 1 1 1
## 744506 0 0 0 0 0 0 1 1 1 0 1 1
## 744507 0 0 0 0 0 0 1 1 1 1 1 0
## 745236 0 0 0 0 0 0 1 1 0 0 0 1
## 745962 0 0 0 1 0 1 0 1 1 0 0 0
## 745963 0 0 0 0 0 0 1 0 1 0 0 0
## 746691 1 0 0 1 1 0 0 0 0 0 0 1
## 746692 0 0 0 1 0 0 1 1 0 0 0 0
## 746693 1 0 0 1 1 1 0 0 0 0 0 0
## 747421 0 0 0 1 0 0 0 0 0 1 1 1
## 747423 1 0 0 1 0 0 0 0 0 1 0 1
## 748150 1 0 0 1 0 0 0 0 1 0 1 1
## 748880 0 0 0 0 1 0 0 0 1 0 1 1
## 748881 0 0 0 1 0 0 0 0 0 0 0 0
## 750260 0 0 0 0 0 0 0 0 0 1 0 0
## 753164 0 0 0 0 0 0 0 0 1 1 0 0
## 753165 0 1 0 0 0 0 0 1 1 1 1 1
## 753900 0 0 0 0 1 0 1 1 1 1 1 1
## 753985 0 0 1 0 0 0 1 0 0 0 1 1
## 753986 0 0 0 0 0 0 0 0 1 1 0 0
## 754607 0 0 0 1 0 1 1 0 1 0 1 1
## 754608 0 0 0 0 1 1 1 1 1 1 1 1
## 754715 0 1 1 0 0 0 1 0 0 0 1 1
## 755444 0 0 0 1 1 0 0 1 1 0 0 0
## 755445 0 0 0 0 0 0 0 1 1 0 1 1
## 755446 0 0 0 0 1 0 0 0 0 0 1 1
## 756173 0 0 0 0 0 1 1 1 1 0 1 1
## 756174 0 0 1 0 1 0 0 1 1 0 1 0
## 756175 0 0 0 0 0 1 0 0 0 0 0 0
## 756903 0 0 0 0 0 0 0 0 0 0 0 1
## 757553 1 1 1 0 1 1 1 1 1 1 1 1
## 757563 0 0 0 0 0 0 1 0 0 0 0 0
## 757634 0 0 0 0 0 0 1 0 0 0 1 1
## 758286 0 0 0 0 0 1 1 1 1 1 1 1
## 758363 0 0 0 0 1 0 1 0 1 0 0 1
## 759023 0 1 0 1 0 1 1 1 1 1 1 1
## 759750 1 0 0 0 0 0 1 1 1 1 1 1
## 759751 1 1 1 0 1 1 1 1 1 1 1 1
## 759753 0 0 0 0 1 1 1 1 1 1 1 1
## 760450 1 0 0 1 0 1 0 0 1 1 1 1
## 760548 0 0 1 1 1 1 1 1 1 0 1 1
## 760549 0 0 0 0 1 1 1 0 1 1 1 1
## 761276 0 0 0 0 0 0 0 1 0 0 1 0
## 761277 0 0 1 0 0 0 0 0 0 1 1 1
## 761278 0 0 0 0 0 1 0 1 0 0 1 1
## 761279 0 0 1 1 0 0 1 1 0 0 1 0
## 762006 1 0 1 1 0 1 1 1 1 1 1 1
## 762007 0 0 1 0 0 0 0 1 1 0 1 1
## 762008 0 0 0 0 1 1 1 1 1 1 1 0
## 762736 0 0 0 1 0 0 0 0 1 1 0 0
## 763424 0 0 0 0 0 0 0 1 1 1 1 0
## 763425 0 0 0 0 1 0 0 1 1 1 0 0
## 764886 0 0 0 0 0 0 0 1 0 0 1 0
## 764919 0 0 0 0 0 0 0 1 0 1 0 1
## 764921 0 0 0 0 0 0 0 1 0 0 0 0
## 765617 0 0 0 0 0 1 1 1 1 1 1 1
## 765649 0 1 0 0 0 0 0 1 1 0 1 1
## 765650 0 0 0 0 0 0 0 0 0 0 1 1
## 766345 0 0 0 0 0 0 0 1 1 1 0 0
## 766346 0 0 0 0 0 0 1 0 0 0 0 1
## 766378 0 0 0 0 1 1 1 1 0 1 1 1
## 766379 0 0 0 0 0 0 1 1 0 0 0 1
## 767073 0 0 0 0 0 0 0 0 0 1 1 0
## 767107 0 0 1 0 0 1 0 1 0 0 0 1
## 768562 1 0 0 0 0 0 1 1 0 0 1 1
## 768564 0 1 1 0 0 0 0 0 0 0 1 0
## 769291 0 0 0 0 0 0 0 0 1 1 1 0
## 769292 1 0 0 1 1 0 1 1 0 0 1 1
## 770017 0 0 0 0 0 0 0 0 0 0 1 1
## 770020 0 0 0 0 1 0 1 0 1 1 0 1
## 770746 0 0 0 0 0 0 1 1 1 1 1 0
## 770747 0 0 0 0 0 0 0 1 1 1 1 1
## 770748 0 0 0 0 0 0 0 0 1 1 1 1
## 771412 0 0 0 0 0 0 0 0 0 0 0 1
## 771474 0 0 0 0 0 0 1 1 1 1 1 1
## 771475 0 1 1 0 0 1 1 1 1 1 1 1
## 771476 0 0 0 0 0 0 0 1 1 1 1 1
## 771477 0 0 0 0 0 0 0 1 0 0 0 0
## 772205 0 1 1 0 0 0 0 1 1 0 1 1
## 772958 0 0 1 1 1 1 0 1 0 1 1 0
## 773687 0 1 0 1 1 1 1 1 1 0 1 1
## 773688 0 0 0 0 0 1 1 1 1 1 0 1
## 774416 0 0 0 0 1 1 0 1 1 0 1 1
## 774417 1 0 0 1 1 1 1 1 1 1 1 1
## 775145 0 1 0 1 1 1 1 1 1 1 1 1
## 775146 1 0 1 1 1 1 1 1 1 1 1 1
## 775147 0 0 0 0 1 0 1 1 1 1 0 1
## 775875 0 0 0 1 0 1 0 0 1 0 1 1
## 775876 0 0 0 0 1 1 1 1 1 1 1 1
## 779475 0 0 0 0 0 0 0 0 0 1 0 0
## 782412 1 1 1 1 1 1 1 1 1 1 1 1
## 782413 1 1 1 1 0 1 1 1 1 1 1 1
## 783150 0 0 1 1 1 1 1 1 1 1 1 1
## 783151 1 1 0 0 1 1 1 1 1 0 1 1
## 783152 0 1 1 1 1 1 1 1 1 1 1 1
## 783872 1 1 1 1 1 1 1 1 1 1 1 1
## 783873 1 1 1 1 1 1 1 1 1 1 1 1
## 783880 0 0 1 0 0 0 0 0 0 0 0 1
## 783881 0 0 0 1 0 0 1 1 0 0 0 0
## 783882 1 1 1 1 1 1 1 1 1 1 1 1
## 783883 0 1 1 1 1 1 1 1 1 1 1 1
## 784603 1 1 1 1 1 1 1 1 1 1 1 1
## 785332 0 1 1 1 1 1 1 1 1 1 1 1
## 789709 0 0 0 1 0 0 0 0 0 0 0 0
## 797711 0 0 0 0 0 0 0 0 0 1 1 1##Filter cells with more than 6 years of sampling
CellsMore6yrs <- CoeflaSimpAnual %>%
group_by(Site.k, Year.t) %>%
summarise(Lists = n()) %>%
group_by(Site.k) %>%
summarise(List_Year = n()) %>%
filter(List_Year >= 6)## `summarise()` has grouped output by 'Site.k'. You can override using the
## `.groups` argument.#Select only cells with more than 6 years of data (from 3736 to 2695)
CoeflaSimpAnual = CoeflaSimpAnual[CoeflaSimpAnual$Site.k %in% CellsMore6yrs$Site.k, ]
#Chronological order
CoeflaSimpAnual <- CoeflaSimpAnual[order(CoeflaSimpAnual$Site.k, CoeflaSimpAnual$Year.t),]
ggplot(CoeflaSimpAnual, aes(x = Year.t, group = Site.k))+
geom_line(aes(y = log(Observed.year.y), color = SpRichness), linetype = "dashed")+
scale_color_gradient2(low = alpha("#4575b4", 0.95),
mid = alpha("#ffffbf", 0.95), midpoint = 73,
high = alpha("#d73027", 0.95))+
geom_point(aes(y = log(Observed.year.y)), size = 0.1, shape = 2)+
scale_y_continuous(expand = c(0,0))+
scale_x_continuous(expand = c(0,0), breaks = c(2012, 2014, 2016, 2018, 2020, 2022))+
labs(x = "Observation date",
y = "Log(mean observed)",
color = "Species richness")+
theme_classic()+
theme(legend.position = "bottom")
onlye some sites to represent the process and run the model
PopID<- unique(CoeflaSimpAnual$Site.k)
#308 cells for conducting the example
#empty list per PopID
Coefla_Parm <- data.frame(cell = character(),
mu.hat = numeric(),
theta.hat = numeric(),
beta.hat = numeric(),
tau.hat = numeric())
for(PopID in unique(CoeflaSimpAnual$Site.k)){
dt_tmp <- CoeflaSimpAnual[CoeflaSimpAnual$Site.k == PopID, ]
Observed.t <- dt_tmp$Observed.year.y
Time.t <- dt_tmp$Year.t
#-------- Log-transform the observations to carry all the calculations in this program--------------------#
log.obs <- log(Observed.t);
#---------------------------------------------------------------------------------------------------------#
#----------------------------------------------------------------------
# PARAMETER ESTIMATION, PARAMETRIC BOOTSTRAP AND PREDICTIONS
#----------------------------------------------------------------------
# Before doing the calculations, the user has to specify ONLY the following 4 options:
# 1. Do you want to compute the ML estimates or the REML estimates?
method <- "REML" # alternatively, set method <- "ML"
# 2. Do you want to plot the predictions?
pred.plot <- "FALSE" # Set it to "FALSE" if you do not want to plot the predictions
# 3. Do you want to plot the parametric bootstrap distribution of the estimates?
pboot.plot <- "FALSE" # Set it to "FALSE" if you do not want to plot the bootstrap distribution of the estimates
# 4. How many bootstrap replicates?
NBoot <- 10; # 10 just to try, for formal results use 1000 or more
#-------------------------------------------------------------------------------------------------#
# 5. RUN THE FOLLOWING LINE OF CODE TO COMPUTE THE ESTIMATES, PREDICTIONS,
# AND RUN A PARAMETRIC BOOTSTRAP. THE USER DOES NOT NEED TO MODIFY THIS LINE OF CODE.
# THE OUTPUT OF THE FUNCTION 'ROUSS.CALCS' IS A LIST AND THE USER CAN RETRIEVE EACH OF THE
# LIST ELEMENTS PRINTED AND SAVED IN THE OBJECT "all.results".
# THE 95% PARAMETRIC BOOTSTRAP FOR BOTH, THE PARAMETERS AND THE PREDICTIONS ARE COMPUTED
# BY THE FUNCTION "ROUSS.CALCS"
#
#-------------------------------------------------------------------------------------------------#
all.results <- ROUSS.CALCS(Yobs=log.obs,Tvec=Time.t, pmethod=method, nboot=NBoot, plot.pred=pred.plot, plot.bootdists = pboot.plot)
Coefla_Parm <- rbind(Coefla_Parm, c(cell = PopID, all.results$parms.est))
print(all.results$parms.est)
}
#It changes the column name, but save the estimates!!
colnames(Coefla_Parm) = c("cell", "mu.hat", "theta.hat", "beta.hat", "tau.hat")
saveRDS(Coefla_Parm, "Coefla_ParametersOUSS_abundance.rds")\(\hat{\mu}\) is the mean stationary log-abundance estimate for each site!
We can also predict the estimated abundance for each year:
PopID<- unique(CoeflaSimpAnual$Site.k)
GSS.Ban <- list()
for(PopID in unique(CoeflaSimpAnual$Site.k)){
ParmPop <- Coefla_Parm %>%
filter(cell == PopID) %>%
dplyr::select(mu.hat, theta.hat, beta.hat, tau.hat)
dt_tmp <- CoeflaSimpAnual[CoeflaSimpAnual$Site.k == PopID, ]
Observed.t <- dt_tmp$Observed.year.y
Tvec <- dt_tmp$Year.t
#-------- Log-transform the observations to carry all the calculations in this program--------------------#
Yobs <- log(Observed.t);
ROUSS.predict(parms = as.numeric(ParmPop), Yobs = Yobs, Tvec = Tvec, plot.it = "TRUE")
PRBan.parms <- c(as.numeric(ParmPop))
#Parameteric Bootstrap function and confidence intervals
B = NBoot
tt <- Tvec-Tvec[1]
long.t <- tt[1]:max(tt)
nparms <- length(PRBan.parms);
preds.boot1<- matrix(0,nrow=B,ncol=length(tt))
preds.boot2<- matrix(0,nrow=B,ncol=length(long.t))
boot.remles <- matrix(0,nrow=B,ncol=nparms+1);
all.sims <- ROUSS.sim(nsims=B,parms=PRBan.parms,Tvec=Tvec);
all.preds <- ROUSS.predict(parms=PRBan.parms, Yobs=Yobs,Tvec=Tvec,plot.it="FALSE")
reml.preds<- all.preds[[1]][,2]
reml.longpreds <- all.preds[[2]][,2]
for(b in 1:B ){
bth.timeseries <- all.sims[,b];
remles.out <- ROUSS.REML(Yobs=bth.timeseries,Tvec=Tvec, parms.guess=PRBan.parms);
boot.remles[b,] <- c(remles.out$remls, remles.out$lnLhat);
all.bootpreds <- ROUSS.predict(parms=remles.out$remls, Yobs=bth.timeseries,Tvec=Tvec,plot.it="FALSE");
preds.boot1[b,] <- all.bootpreds[[1]][,2]
preds.boot2[b,] <- all.bootpreds[[2]][,2]
}
CIs.mat <- apply(boot.remles,2,FUN=function(x){quantile(x,probs=c(0.025,0.975))});
CIs.mat <- rbind(CIs.mat[1,1:4],PRBan.parms,CIs.mat[2,1:4]);
rownames(CIs.mat) <- c("2.5%","REMLE","97.5%");
colnames(CIs.mat) <- c("mu", "theta","betasq","tausq");
preds.CIs1 <- apply(preds.boot1,2,FUN=function(x){quantile(x,probs=c(0.025,0.975))});
mean.boots <- apply(preds.boot1,2,FUN=function(x){quantile(x,probs=0.50)})
preds.CIs1 <- t(rbind(Tvec,reml.preds-(mean.boots-preds.CIs1[1,]), reml.preds, reml.preds+(preds.CIs1[2,]-mean.boots)));
colnames(preds.CIs1) <- c("Year","2.5%","REMLE","97.5%");
preds.CIs2 <- apply(preds.boot2,2,FUN=function(x){quantile(x,probs=c(0.025,0.975))});
mean.boots2 <- apply(preds.boot2,2,FUN=function(x){quantile(x,probs=0.50)})
preds.CIs2 <- t(rbind(long.t+Tvec[1],reml.longpreds-(mean.boots2-preds.CIs2[1,]), reml.longpreds, reml.longpreds+(preds.CIs2[2,]-mean.boots2)));
#preds.CIs2 <- apply(preds.boot2,2,FUN=function(x){quantile(x,probs=c(0.025,0.5,0.975))});
#preds.CIs2 <- t(rbind(long.t+Tvec[1],preds.CIs2));
colnames(preds.CIs2) <- c("Year","2.5%","REMLE","97.5%");
#pred.CIs2 <- cbind(preds.CIs2[,1], reml.longpreds-(preds.CIs2[,3]-preds.CIs2[,2]),reml.longpreds,reml.longpreds+(preds.CIs2[,4]-preds.CIs2[,3]))
boot.list <- list(boot.remles = boot.remles, CIs.mat = CIs.mat, preds.CIs1 = preds.CIs1,
preds.CIs2=preds.CIs2)
boot.list$preds.CIs1
df.boot <- as.data.frame(cbind(PopID, boot.list$preds.CIs1))
GSS.Ban <- c(GSS.Ban, list(df.boot))
}
#extract estimates
GSS.BanAll <- bind_rows(GSS.Ban)
colnames(GSS.BanAll) <-c("Site.k", "Year.t","Low","REMLE","High");
GSS.BanAll <- GSS.BanAll %>%
mutate_at(c("Year.t","Low","REMLE","High"), as.numeric) %>%
left_join(CoeflaSimpAnual, by = c("Site.k", "Year.t"))
saveRDS(GSS.BanAll, "Coefla_predictedAbundance_REMLE_OUSS.rds")GSS.BanAll <- readRDS("Coefla_predictedAbundance_REMLE_OUSS.rds")
ggplot(GSS.BanAll, aes(x = Year.t))+
geom_line(aes(y = log(Observed.year.y), group = Site.k, color = SpRichness), linetype = "dotted")+
geom_point(aes(y = log(Observed.year.y), fill = SpRichness), shape = 21)+
geom_line(aes(y = log(REMLE), group = Site.k, color = SpRichness), linetype = "solid")+
geom_point(aes(y = log(REMLE), fill = SpRichness), shape = 22)+
scale_color_gradient2(low = alpha("#4575b4", 0.95),
mid = alpha("#ffffbf", 0.95), midpoint = 73,
high = alpha("#d73027", 0.95))+
scale_fill_gradient2(low = alpha("#4575b4", 0.95),
mid = alpha("#ffffbf", 0.95), midpoint = 73,
high = alpha("#d73027", 0.95))+
labs(x = "t",
y = "log(Observed / estimated abundance)",
color = "Species richness",
fill = "Species richness",
title = "62 sites (hexagonal cells)")+
theme_classic()+
theme(legend.position = c(0.8,0.8),
legend.direction = "horizontal")## Warning: A numeric `legend.position` argument in `theme()` was deprecated in ggplot2
## 3.5.0.
## ℹ Please use the `legend.position.inside` argument of `theme()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
ggplot(GSS.BanAll, aes(x = SpRichness, fill = SpRichness))+
geom_point(aes(y = REMLE), shape = 22)+
geom_point(aes(y = Observed.year.y), shape = 21)+
geom_smooth(aes(x = SpRichness, y = REMLE), method = "lm", linetype = "solid")+
geom_smooth(aes(x = SpRichness, y = Observed.year.y), method = "lm", linetype = "dashed")+
scale_fill_gradient2(low = alpha("#4575b4", 0.95),
mid = alpha("#ffffbf", 0.95), midpoint = 73,
high = alpha("#d73027", 0.95))+
labs(x = "Species richness",
y = "Observed / estimated abundance",
fill = "Species richness",
title = "Relationship of abundance and species richness")+
theme_classic()+
theme(legend.position = c(0.75,0.75),
legend.direction = "horizontal")## `geom_smooth()` using formula = 'y ~ x'## Warning: The following aesthetics were dropped during statistical transformation: fill.
## ℹ This can happen when ggplot fails to infer the correct grouping structure in
## the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
## variable into a factor?## `geom_smooth()` using formula = 'y ~ x'## Warning: The following aesthetics were dropped during statistical transformation: fill.
## ℹ This can happen when ggplot fails to infer the correct grouping structure in
## the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
## variable into a factor?
We can map our data
GSS.map.df <- GSS.BanAll %>%
group_by(Site.k, Year.t, REMLE) %>%
summarise(SpRichness = mean(SpRichness),
lon = mean(longitude),
lat = mean(latitude))## `summarise()` has grouped output by 'Site.k', 'Year.t'. You can override using
## the `.groups` argument.ggplot() +
geom_sf(data = world1)+
geom_point(data = GSS.map.df,
aes(x = lon, y = lat, size = REMLE, fill = SpRichness), alpha = 0.5, shape = 21)+
scale_fill_gradient2(low = alpha("#4575b4", 0.75),
mid = alpha("#ffffbf", 0.75), midpoint = 73,
high = alpha("#d73027", 0.75))+
coord_sf(xlim = c(-91, -59.4),
ylim = c(8, 27)) +
scale_x_continuous(breaks = seq(from = -90, to = -60, by = 15))+
labs(y = "Latitude",
x = "Longitude",
tag = "",
title = "Caribbean (high eBird coverage-based completeness)",
subtitle = "Abundance estimation using GSS(OUSS) models for Bananaquit",
fill = "Species richness",
size = "Stationary abundance estimation") +
theme_classic()+
theme(legend.position = "bottom",
legend.direction = "horizontal")+
geom_magnify(data = GSS.map.df, aes(x = lon, y = lat),
from = c(xmin = -86.2, xmax = -81.2, ymin = 8, ymax = 13.5),
to = c(xmin = -75, xmax = -62, ymin = 10, ymax = 25),
shadow = TRUE) ## Warning in layer_sf(geom = ggproto(NULL, GeomMagnify), mapping = mapping, :
## Ignoring unknown aesthetics: x and y## Warning in .setMask(mask$f, mask$ref): Ignored alpha mask (not supported on
## this device)
## Warning in .setMask(mask$f, mask$ref): Ignored alpha mask (not supported on
## this device)
## Warning in .setMask(mask$f, mask$ref): Ignored alpha mask (not supported on
## this device)
## Warning in .setMask(mask$f, mask$ref): Ignored alpha mask (not supported on
## this device)
## Warning in .setMask(mask$f, mask$ref): Ignored alpha mask (not supported on
## this device)
## Warning in .setMask(mask$f, mask$ref): Ignored alpha mask (not supported on
## this device)
## Warning in .setMask(mask$f, mask$ref): Ignored alpha mask (not supported on
## this device)
## Warning in .setMask(mask$f, mask$ref): Ignored alpha mask (not supported on
## this device)
## Warning in .setMask(mask$f, mask$ref): Ignored alpha mask (not supported on
## this device)
## Warning in .setMask(mask$f, mask$ref): Ignored alpha mask (not supported on
## this device)
## Warning in .setMask(mask$f, mask$ref): Ignored alpha mask (not supported on
## this device)
## Warning in .setMask(mask$f, mask$ref): Ignored alpha mask (not supported on
## this device)
## Warning in .setMask(mask$f, mask$ref): Ignored alpha mask (not supported on
## this device)
## Warning in .setMask(mask$f, mask$ref): Ignored alpha mask (not supported on
## this device)
## Warning in .setMask(mask$f, mask$ref): Ignored alpha mask (not supported on
## this device)
## Warning in .setMask(mask$f, mask$ref): Ignored alpha mask (not supported on
## this device)
## Warning in .setMask(mask$f, mask$ref): Ignored alpha mask (not supported on
## this device)
## Warning in .setMask(mask$f, mask$ref): Ignored alpha mask (not supported on
## this device)
## Warning in .setMask(mask$f, mask$ref): Ignored alpha mask (not supported on
## this device)
## Warning in .setMask(mask$f, mask$ref): Ignored alpha mask (not supported on
## this device)
## Warning in .setMask(mask$f, mask$ref): Ignored alpha mask (not supported on
## this device)
## Warning in .setMask(mask$f, mask$ref): Ignored alpha mask (not supported on
## this device)
## Warning in .setMask(mask$f, mask$ref): Ignored alpha mask (not supported on
## this device)
## Warning in .setMask(mask$f, mask$ref): Ignored alpha mask (not supported on
## this device)
## Warning in .setMask(mask$f, mask$ref): Ignored alpha mask (not supported on
## this device)
## Warning in .setMask(mask$f, mask$ref): Ignored alpha mask (not supported on
## this device)
## Warning in .setMask(mask$f, mask$ref): Ignored alpha mask (not supported on
## this device)
## Warning in .setMask(mask$f, mask$ref): Ignored alpha mask (not supported on
## this device)
End of this document. Orlando.
A previous approach (reducing the observer error parameter and applying Data Cloning)
We tried with the Poisson process (adjust to the length of the time series). Here is the GSS.dc model:
StochGSS2.dc <- function(){
# Priors on model parameters. Priors are DC1 in Lele et al (2007)
a1 ~ dnorm(0,1); # constant, the population growth rate.
c1 ~ dunif(-1,1); # constant, the density dependence parameter.
sig1 ~ dlnorm(-0.5,10); #variance parameter of stochastic environment (process noise) in the system
stovar1 <- 1/pow(sig1,2)
for(k in 1:K){
# Simulate trajectory that depends on the previous
mean_X1[1,k] <- a1/(1-c1) # Expected value of the first realization of the process
# this is drawn from the stationary distribution of the process
# Equation 14 (main text) and A.4 in Appendix of Dennis et al 2006
Varno1[k] <- pow(sig1,2)/(1-pow(c1,2)) #. Equation A.5 in Appendix of Dennis et al 2006
# Updating the state: Stochastic process for 29 steps
X1[1,k]~dnorm(mean_X1[1,k], 1/Varno1[k]); #first estimation of population
#iteration of the GSS model in the data
for (i in 2:12) {
mean_X1[i,k] <- a1 + c1 * X1[(i - 1),k]
X1[i,k] ~ dnorm(mean_X1[i,k], stovar1) # Et is included here since a+cX(t-1) + Et ~ Normal(a+cX(t-1),sigma^2)
}
# Updating the observations, from the counts
for (i in 1:12) {
Y1[i,k] ~ dpois(exp(X1[i,k])) #Here we can plug the probability of detection, p*exp(X1[i,k]), and we could also estimate it from the data!
}
}
}Now, we select the data, first Costa Rica - Alejuela (CR-A), which has complete data 2012-2023.
table(Pop$state_code, Pop$year)
#Costa Rica - Alejuela
CR.A <- Pop %>%
dplyr::filter(state_code == "CR-A") %>%
dplyr::select(mu_count)
CR <- round(CR.A$mu_count, 0)
datalistGSS.CR.dc <- list(K = 1,
Y1 = dcdim(data.matrix(CR))) # Not need to be in LOG!!!
dcrun.GSS2.CR <- dc.fit(data = datalistGSS.CR.dc,
params = c("a1", "c1", "sig1"), #extract from the model
model = StochGSS2.dc,
n.clones = c(5,10,50,100,150),
multiply = "K",
n.chains = 10,
n.adapt = 10000,
n.update = 10000,
thin = 10,
n.iter = 100000)
summary(dcrun.GSS2.CR);
dcdiag(dcrun.GSS2.CR) Now, Aruba (Lesser Antilles), missing the two first years of data, complete data 2014-2023
A <- Pop %>%
dplyr::filter(state_code == "AW-") %>%
dplyr::select(mu_count)
AW <- c(round(A$mu_count, 0)) #the
#and have to adjust the model
StochGSS2.AW.dc <- function(){
# Priors on model parameters. Priors are DC1 in Lele et al (2007)
a1 ~ dnorm(0,1); # constant, the population growth rate.
c1 ~ dunif(-1,1); # constant, the density dependence parameter.
sig1 ~ dlnorm(-0.5,10); #variance parameter of stochastic environment (process noise) in the system
stovar1 <- 1/pow(sig1,2)
for(k in 1:K){
# Simulate trajectory that depends on the previous
mean_X1[1,k] <- a1/(1-c1) # Expected value of the first realization of the process
# this is drawn from the stationary distribution of the process
# Equation 14 (main text) and A.4 in Appendix of Dennis et al 2006
Varno1[k] <- pow(sig1,2)/(1-pow(c1,2)) #. Equation A.5 in Appendix of Dennis et al 2006
# Updating the state: Stochastic process for 29 steps
X1[1,k]~dnorm(mean_X1[1,k], 1/Varno1[k]); #first estimation of population
#iteration of the GSS model in the data
for (i in 2:10) {
mean_X1[i,k] <- a1 + c1 * X1[(i - 1),k]
X1[i,k] ~ dnorm(mean_X1[i,k], stovar1) # Et is included here since a+cX(t-1) + Et ~ Normal(a+cX(t-1),sigma^2)
}
# Updating the observations, from the counts
for (i in 1:10) {
Y1[i,k] ~ dpois(exp(X1[i,k])) #Here we can plug the probability of detection, p*exp(X1[i,k]), and we could also estimate it from the data!
}
}
}
datalistGSS.AW.dc <- list(K = 1,
Y1 = dcdim(data.matrix(AW))) # Not need to be in LOG!!!
dcrun.GSS2.AW <- dc.fit(data = datalistGSS.AW.dc,
params = c("a1", "c1", "sig1"), #extract from the model
model = StochGSS2.AW.dc,
n.clones = c(5,10,50,100,150),
multiply = "K",
n.chains = 10,
n.adapt = 10000,
n.update = 10000,
thin = 10,
n.iter = 100000)
summary(dcrun.GSS2.AW);
dcdiag(dcrun.GSS2.AW) Let’s try Puerto Rico - San Juan, which has NA in 2013 and 2014. I will run from 2015 to 2023
PR.SJ <- Pop %>%
dplyr::filter(state_code == "PR-SJ") %>%
dplyr::select(mu_count)
PR <- c(round(PR.SJ$mu_count, 0)) #the
#and have to adjust the model
StochGSS2.PR.dc <- function(){
# Priors on model parameters. Priors are DC1 in Lele et al (2007)
a1 ~ dnorm(0,1); # constant, the population growth rate.
c1 ~ dunif(-1,1); # constant, the density dependence parameter.
sig1 ~ dlnorm(-0.5,10); #variance parameter of stochastic environment (process noise) in the system
stovar1 <- 1/pow(sig1,2)
for(k in 1:K){
# Simulate trajectory that depends on the previous
mean_X1[1,k] <- a1/(1-c1) # Expected value of the first realization of the process
# this is drawn from the stationary distribution of the process
# Equation 14 (main text) and A.4 in Appendix of Dennis et al 2006
Varno1[k] <- pow(sig1,2)/(1-pow(c1,2)) #. Equation A.5 in Appendix of Dennis et al 2006
# Updating the state: Stochastic process for 29 steps
X1[1,k]~dnorm(mean_X1[1,k], 1/Varno1[k]); #first estimation of population
#iteration of the GSS model in the data
for (i in 2:10) {
mean_X1[i,k] <- a1 + c1 * X1[(i - 1),k]
X1[i,k] ~ dnorm(mean_X1[i,k], stovar1) # Et is included here since a+cX(t-1) + Et ~ Normal(a+cX(t-1),sigma^2)
}
# Updating the observations, from the counts
for (i in 1:10) {
Y1[i,k] ~ dpois(exp(X1[i,k])) #Here we can plug the probability of detection, p*exp(X1[i,k]), and we could also estimate it from the data!
}
}
}
datalistGSS.PR.dc <- list(K = 1,
Y1 = dcdim(data.matrix(PR))) # Not need to be in LOG!!!
dcrun.GSS2.PR <- dc.fit(data = datalistGSS.PR.dc,
params = c("a1", "c1", "sig1"), #extract from the model
model = StochGSS2.PR.dc,
n.clones = c(5,10,50,100,150),
multiply = "K",
n.chains = 10,
n.adapt = 10000,
n.update = 10000,
thin = 10,
n.iter = 100000)
summary(dcrun.GSS2.PR);
dcdiag(dcrun.GSS2.PR) Let’s try Venezuela - Miranda, which has data from 2016 to 2023
VE.M <- Pop %>%
dplyr::filter(state_code == "VE-M") %>%
dplyr::select(mu_count)
VE <- c(round(VE.M$mu_count, 0)) #the
#and have to adjust the model
StochGSS2.VE.dc <- function(){
# Priors on model parameters. Priors are DC1 in Lele et al (2007)
a1 ~ dnorm(0,1); # constant, the population growth rate.
c1 ~ dunif(-1,1); # constant, the density dependence parameter.
sig1 ~ dlnorm(-0.5,10); #variance parameter of stochastic environment (process noise) in the system
stovar1 <- 1/pow(sig1,2)
for(k in 1:K){
# Simulate trajectory that depends on the previous
mean_X1[1,k] <- a1/(1-c1) # Expected value of the first realization of the process
# this is drawn from the stationary distribution of the process
# Equation 14 (main text) and A.4 in Appendix of Dennis et al 2006
Varno1[k] <- pow(sig1,2)/(1-pow(c1,2)) #. Equation A.5 in Appendix of Dennis et al 2006
# Updating the state: Stochastic process for 29 steps
X1[1,k]~dnorm(mean_X1[1,k], 1/Varno1[k]); #first estimation of population
#iteration of the GSS model in the data
for (i in 2:9) {
mean_X1[i,k] <- a1 + c1 * X1[(i - 1),k]
X1[i,k] ~ dnorm(mean_X1[i,k], stovar1) # Et is included here since a+cX(t-1) + Et ~ Normal(a+cX(t-1),sigma^2)
}
# Updating the observations, from the counts
for (i in 1:9) {
Y1[i,k] ~ dpois(exp(X1[i,k])) #Here we can plug the probability of detection, p*exp(X1[i,k]), and we could also estimate it from the data!
}
}
}
datalistGSS.VE.dc <- list(K = 1,
Y1 = dcdim(data.matrix(VE))) # Not need to be in LOG!!!
dcrun.GSS2.VE <- dc.fit(data = datalistGSS.VE.dc,
params = c("a1", "c1", "sig1"), #extract from the model
model = StochGSS2.VE.dc,
n.clones = c(5,10,50,100,150),
multiply = "K",
n.chains = 10,
n.adapt = 10000,
n.update = 10000,
thin = 10,
n.iter = 100000)
summary(dcrun.GSS2.VE);
dcdiag(dcrun.GSS2.VE)Let’s try Grenada - Saint Andrew, which has six years of data
GR.SA <- Pop %>%
dplyr::filter(state_code == "GD-01") %>%
dplyr::select(mu_count)
GR <- c(round(GR.SA$mu_count, 0)) #the
#and have to adjust the model
StochGSS2.GR.dc <- function(){
# Priors on model parameters. Priors are DC1 in Lele et al (2007)
a1 ~ dnorm(0,1); # constant, the population growth rate.
c1 ~ dunif(-1,1); # constant, the density dependence parameter.
sig1 ~ dlnorm(-0.5,10); #variance parameter of stochastic environment (process noise) in the system
stovar1 <- 1/pow(sig1,2)
for(k in 1:K){
# Simulate trajectory that depends on the previous
mean_X1[1,k] <- a1/(1-c1) # Expected value of the first realization of the process
# this is drawn from the stationary distribution of the process
# Equation 14 (main text) and A.4 in Appendix of Dennis et al 2006
Varno1[k] <- pow(sig1,2)/(1-pow(c1,2)) #. Equation A.5 in Appendix of Dennis et al 2006
# Updating the state: Stochastic process for 29 steps
X1[1,k]~dnorm(mean_X1[1,k], 1/Varno1[k]); #first estimation of population
#iteration of the GSS model in the data
for (i in 2:6) {
mean_X1[i,k] <- a1 + c1 * X1[(i - 1),k]
X1[i,k] ~ dnorm(mean_X1[i,k], stovar1) # Et is included here since a+cX(t-1) + Et ~ Normal(a+cX(t-1),sigma^2)
}
# Updating the observations, from the counts
for (i in 1:6) {
Y1[i,k] ~ dpois(exp(X1[i,k])) #Here we can plug the probability of detection, p*exp(X1[i,k]), and we could also estimate it from the data!
}
}
}
datalistGSS.GR.dc <- list(K = 1,
Y1 = dcdim(data.matrix(GR))) # Not need to be in LOG!!!
dcrun.GSS2.GR <- dc.fit(data = datalistGSS.GR.dc,
params = c("a1", "c1", "sig1"), #extract from the model
model = StochGSS2.GR.dc,
n.clones = c(5,10,50,100,150),
multiply = "K",
n.chains = 10,
n.adapt = 10000,
n.update = 10000,
thin = 10,
n.iter = 100000)
summary(dcrun.GSS2.GR);
dcdiag(dcrun.GSS2.GR)And Jamaica - Saint Andrew, which has four years of data
JM.SA <- Pop %>%
dplyr::filter(state_code == "JM-02") %>%
dplyr::select(mu_count)
JM <- c(round(JM.SA$mu_count, 0)) #the
#and have to adjust the model
StochGSS2.JM.dc <- function(){
# Priors on model parameters. Priors are DC1 in Lele et al (2007)
a1 ~ dnorm(0,1); # constant, the population growth rate.
c1 ~ dunif(-1,1); # constant, the density dependence parameter.
sig1 ~ dlnorm(-0.5,10); #variance parameter of stochastic environment (process noise) in the system
stovar1 <- 1/pow(sig1,2)
for(k in 1:K){
# Simulate trajectory that depends on the previous
mean_X1[1,k] <- a1/(1-c1) # Expected value of the first realization of the process
# this is drawn from the stationary distribution of the process
# Equation 14 (main text) and A.4 in Appendix of Dennis et al 2006
Varno1[k] <- pow(sig1,2)/(1-pow(c1,2)) #. Equation A.5 in Appendix of Dennis et al 2006
# Updating the state: Stochastic process for 29 steps
X1[1,k]~dnorm(mean_X1[1,k], 1/Varno1[k]); #first estimation of population
#iteration of the GSS model in the data
for (i in 2:4) {
mean_X1[i,k] <- a1 + c1 * X1[(i - 1),k]
X1[i,k] ~ dnorm(mean_X1[i,k], stovar1) # Et is included here since a+cX(t-1) + Et ~ Normal(a+cX(t-1),sigma^2)
}
# Updating the observations, from the counts
for (i in 1:4) {
Y1[i,k] ~ dpois(exp(X1[i,k])) #Here we can plug the probability of detection, p*exp(X1[i,k]), and we could also estimate it from the data!
}
}
}
datalistGSS.JM.dc <- list(K = 1,
Y1 = dcdim(data.matrix(JM))) # Not need to be in LOG!!!
dcrun.GSS2.JM <- dc.fit(data = datalistGSS.JM.dc,
params = c("a1", "c1", "sig1"), #extract from the model
model = StochGSS2.JM.dc,
n.clones = c(5,10,50,100,150),
multiply = "K",
n.chains = 10,
n.adapt = 10000,
n.update = 10000,
thin = 10,
n.iter = 100000)
summary(dcrun.GSS2.JM);
dcdiag(dcrun.GSS2.JM)Now the null:
H0 <- round(Pop$mu_count, 0)
#and have to adjust the model
StochGSS2.H0.dc <- function(){
# Priors on model parameters. Priors are DC1 in Lele et al (2007)
a1 ~ dnorm(0,1); # constant, the population growth rate.
c1 ~ dunif(-1,1); # constant, the density dependence parameter.
sig1 ~ dlnorm(-0.5,10); #variance parameter of stochastic environment (process noise) in the system
stovar1 <- 1/pow(sig1,2)
for(k in 1:K){
# Simulate trajectory that depends on the previous
mean_X1[1,k] <- a1/(1-c1) # Expected value of the first realization of the process
# this is drawn from the stationary distribution of the process
# Equation 14 (main text) and A.4 in Appendix of Dennis et al 2006
Varno1[k] <- pow(sig1,2)/(1-pow(c1,2)) #. Equation A.5 in Appendix of Dennis et al 2006
# Updating the state: Stochastic process for 29 steps
X1[1,k]~dnorm(mean_X1[1,k], 1/Varno1[k]); #first estimation of population
#iteration of the GSS model in the data
for (i in 2:51) {
mean_X1[i,k] <- a1 + c1 * X1[(i - 1),k]
X1[i,k] ~ dnorm(mean_X1[i,k], stovar1) # Et is included here since a+cX(t-1) + Et ~ Normal(a+cX(t-1),sigma^2)
}
# Updating the observations, from the counts
for (i in 1:51) {
Y1[i,k] ~ dpois(exp(X1[i,k])) #Here we can plug the probability of detection, p*exp(X1[i,k]), and we could also estimate it from the data!
}
}
}
datalistGSS.H0.dc <- list(K = 1,
Y1 = dcdim(data.matrix(H0))) # Not need to be in LOG!!!
dcrun.GSS2.H0 <- dc.fit(data = datalistGSS.H0.dc,
params = c("a1", "c1", "sig1"), #extract from the model
model = StochGSS2.H0.dc,
n.clones = c(5,10,50,100,150),
multiply = "K",
n.chains = 10,
n.adapt = 10000,
n.update = 10000,
thin = 10,
n.iter = 100000)
summary(dcrun.GSS2.H0);
dcdiag(dcrun.GSS2.H0)Let’s combine the results
MLE <- as.data.frame(rbind(coef(dcrun.GSS2.H0),
coef(dcrun.GSS2.CR),
coef(dcrun.GSS2.VE),
coef(dcrun.GSS2.AW),
coef(dcrun.GSS2.GR),
coef(dcrun.GSS2.PR),
coef(dcrun.GSS2.JM)))
MLE$Pop <- c("Null",
"CR-A",
"VE-M",
"AW-",
"GD-01",
"PR-SJ",
"JM-02")
MLE$N_equ <- exp(-(MLE$a1/(MLE$c1-1)))
print(MLE)
#Some figures
#how many records per site?
Pop %>%
group_by(state_code) %>%
summarise(AllRecords = sum(n_records))Mainland sites
par(mfrow = c(1,2));#create a 2x2 plotting matrix
hist(rnorm(10000, mean = MLE[2,5], sd = MLE[2,3]/(1-MLE[2,2])),
xlab = "Equilibrium density",
main = "Histogram of Costa Rica - Alejuela",
xlim = c(0,10));
abline(v=MLE[2,5], col="blue");
abline(v=mean(CR), col = "red");
legend(5, 1500, legend=c("Abundance (estimated)", "Observation mean"),
col=c("blue", "red"), lty=1, cex=0.8);
text(x = 8, y = 500, "n = 249")
hist(rnorm(10000, mean = MLE[3,5], sd = MLE[3,3]/(1-MLE[3,2])),
xlab = "Equilibrium density",
main = "Histogram of Venezuela - Miranda",
xlim = c(0,10));
abline(v=MLE[3,5], col="blue");
abline(v=mean(VE), col = "red");
text(x = 8, y = 500, "n = 521")Lesser Antilles sites
par(mfrow = c(1,2));#create a 2x2 plotting matrix
hist(rnorm(10000, mean = MLE[4,5], sd = MLE[4,3]/(1-MLE[4,2])),
xlab = "Equilibrium density",
main = "Histogram of Aruba",
xlim = c(0,10));
abline(v=MLE[4,5], col="blue");
abline(v=mean(AW), col = "red");
text(x = 8, y = 500, "n = 229")
hist(rnorm(10000, mean = MLE[5,5], sd = MLE[5,3]/(1-MLE[5,2])),
xlab = "Equilibrium density",
main = "Histogram of Grenada - Saint Andrew",
xlim = c(0,10));
abline(v=MLE[5,5], col="blue");
abline(v=mean(GR), col = "red");
text(x = 8, y = 500, "n = 117")Greater Antilles sites
par(mfrow = c(1,2));#create a 2x2 plotting matrix
hist(rnorm(10000, mean = MLE[6,5], sd = MLE[6,3]/(1-MLE[6,2])),
xlab = "Equilibrium density",
main = "Histogram of Puerto Rico - San Juan",
xlim = c(0,10));
abline(v=MLE[6,5], col="blue");
abline(v=mean(PR), col = "red");
text(x = 8, y = 500, "n = 105")
hist(rnorm(10000, mean = MLE[7,5], sd = MLE[7,3]/(1-MLE[7,2])),
xlab = "Equilibrium density",
main = "Histogram of Jamaica - Saint Andrew",
xlim = c(0,10));
abline(v=MLE[7,5], col="blue");
abline(v=mean(JM), col = "red");
text(x = 8, y = 500, "n = 18")Null hypothesis and relationship
par(mfrow = c(1,2));#create a 2x2 plotting matrix
hist(rnorm(10000, mean = MLE[1,5], sd = MLE[1,3]/(1-MLE[1,2])),
xlab = "Equilibrium density",
main = "Histogram of all the Meta-archipelago of the Caribbean",
xlim = c(0,10));
abline(v=MLE[1,5], col="blue");
abline(v=mean(H0), col = "red");
text(x = 8, y = 500, "n = 1,239")
MLE$S <- c(NA, 149, 76, 67, 51, 39, 31)
plot(x = MLE[c(2:7),6], y = MLE[c(2:7),5], pch = 19, cex = 4,
xlim = c(0,150), ylim = c(0,5),
xlab = "Species richness",
ylab = "Abundance (estimated)",
main = "Hypothetical relationship",
col = c("#F36F20", "#DD692E", "#864268", "#704C76", "#434174","#01458A"));
abline(lm(MLE$N_equ ~ MLE$S), col = "darkgray");
text(x = 25, y = 1, expression('R' ^ 2~"= 0.189; p = 0.215"))Likelihood Ratio test
The parameters for the estimation of the unknown, mean number of Bananaquit per list of ≤2 km and 10-120 minutes in the six sites \(s\) (\(s = CR, VE, AR, GD, PR, JM\)) are \(\mathbf{\Theta} = [a, c, \sigma]\). Several possible scenarios regarding the comparison of these quantities across sites comes to mind. The simple scenario is that all parameters are equal among the meta-archipelago and representing identical parameters \(\mathbf{\Theta}\), thus our working null hypothesis is
\[ {\rm H_{0}} : \mathbf{\Theta}_{CR} = \mathbf{\Theta}_{VE} = \mathbf{\Theta}_{AW} = \mathbf{\Theta}_{GD} = \mathbf{\Theta}_{PR} = \mathbf{\Theta}_{JM} = \mathbf{\Theta} \]
A very contrasting alternative hypothesis (Model1) represent differences among the sites \(s\) in the meta-archipelago, and is \[ {\rm H_{1}} : \mathbf{\Theta}_{CR} \neq \mathbf{\Theta}_{VE} \neq \mathbf{\Theta}_{AW} \neq \mathbf{\Theta}_{GD} \neq \mathbf{\Theta}_{PR} \neq \mathbf{\Theta}_{JM} \]
A third hypothesis is that sub-regions within the meta-archipelago (Mainland vs Lesser Antilles vs Greater Antilles) have similar parameters \(\mathbf{\Theta}_{Subregion}\), and they differ from the others…
\[ {\rm H_{2}} : \mathbf{\Theta}_{CR} = \mathbf{\Theta}_{VE} \neq \mathbf{\Theta}_{AW} = \mathbf{\Theta}_{GD} \neq \mathbf{\Theta}_{PR} = \mathbf{\Theta}_{JM} \] or \[ {\rm H_{2}} : \mathbf{\Theta}_{Mainland} \neq \mathbf{\Theta}_{Lesser Antilles} \neq \mathbf{\Theta}_{Greater Antilles} \]
And another hypothesis is that representing “the island syndrome” (Mainland vs Islands) of difference in parameters \(\mathbf{\Theta}_{Mainland-Island}\):
\[ {\rm H_{3}} : \mathbf{\Theta}_{CR} = \mathbf{\Theta}_{VE} \neq \mathbf{\Theta}_{AW} = \mathbf{\Theta}_{GD} = \mathbf{\Theta}_{PR} = \mathbf{\Theta}_{JM} \] or \[ {\rm H_{3}} : \mathbf{\Theta}_{Mainland} \neq \mathbf{\Theta}_{Islands} \]
The LR test statistics are \[ {\Lambda_{0-1}} = \frac{L_0(\mathbf{\hat\Theta})}{L_1(\mathbf{\hat\Theta}_{CR},\mathbf{\hat\Theta}_{VE},\mathbf{\hat\Theta}_{AW},\mathbf{\hat\Theta}_{GD},\mathbf{\hat\Theta}_{PR},\mathbf{\hat\Theta}_{JM})} \]
\[ {\Lambda_{0-2}} = \frac{L_0(\hat\Theta)}{L_2(\mathbf{\hat\Theta}_{Mainland},\mathbf{\hat\Theta}_{Lesser Antilles},\mathbf{\hat\Theta}_{Greater Antilles})} \]
\[ {\Lambda_{0-3}} = \frac{L_0(\hat\Theta)}{L_3(\mathbf{\hat\Theta}_{Mainland},\mathbf{\hat\Theta}_{Islands})} \]
\[ {\Lambda_{2-1}} = \frac{L_2(\mathbf{\hat\Theta}_{Mainland},\mathbf{\hat\Theta}_{Lesser Antilles},\mathbf{\hat\Theta}_{Greater Antilles})}{L_1(\mathbf{\hat\Theta}_{CR},\mathbf{\hat\Theta}_{VE},\mathbf{\hat\Theta}_{AW},\mathbf{\hat\Theta}_{GD},\mathbf{\hat\Theta}_{PR},\mathbf{\hat\Theta}_{JM})} \]